#include "cuda_runtime.h"
#include "device_launch_parameters.h"

#include <stdio.h>
#include <fstream>
#include <time.h>
#include <cmath>
#include <string>
#include <sm_35_atomic_functions.h>

using namespace std;

#define PI 3.1415926535897932384626433832795
#define RunTime1 15e5
#define RunTime2 2e6
#define Level 57
#define NoC 56//(level-1)
#define NoM 112//2*(level-1)
#define NoC1 48//(level-1)
#define NoM1 96//2*(level-1)
#define freq 50.0
#define M 6
#define N 7
#define Ma 14
#define Na 15
#define Mb 2
#define Nb 3

#define NoSTN 2
#define NLine 1


__device__ void solver(double AMatrix[M][N], double I[M])
{
	double temp[M] = { 0 };
	int k = 0;

	while (k < M - 1)
	{
		temp[k] = AMatrix[k][k];

		for (int j = k; j < N; j++)
		{
			AMatrix[k][j] = AMatrix[k][j] / temp[k];
		}
		for (int i = k + 1; i < M; i++)
		{
			//for (int j = k + 1; j < N; j++)
			//{
			//	AMatrix[i][j] = AMatrix[i][j] - AMatrix[i][k] * AMatrix[k][j];
			//}
			//AMatrix[i][k] = 0;

			if (AMatrix[i][k] != 0)
			{
				for (int j = k + 1; j < N; j++)
				{
					AMatrix[i][j] = AMatrix[i][j] - AMatrix[i][k] * AMatrix[k][j];
				}
			}
		}
		k = k + 1;
	}

	double Adiag = AMatrix[M - 1][M - 1];
	for (int j = M - 1; j <= N - 1; j++)
	{
		AMatrix[M - 1][j] = AMatrix[M - 1][j] / Adiag;
	}

	int kk = M - 1;
	while (kk > 0)
	{
		for (int i = 0; i <= kk - 1; i++)//or for (int i=0;i<kk-1;i++)?
		{
			temp[i] = AMatrix[i][kk];
			for (int j = kk; j <= N - 1; j++)
			{
				AMatrix[i][j] = AMatrix[i][j] - temp[i] * AMatrix[kk][j];
			}
		}
		kk = kk - 1;
	}


	for (int i = 0; i < M; i++)
	{
		I[i] = AMatrix[i][N - 1];
	}
}

__device__ void solverSM(double AMatrix[Ma][Na], double I[Ma])
{
	double temp[Ma] = { 0 };
	int k = 0;

	while (k < Ma - 1)
	{
		temp[k] = AMatrix[k][k];

		for (int j = k; j < Na; j++)
		{
			AMatrix[k][j] = AMatrix[k][j] / temp[k];
		}
		for (int i = k + 1; i < Ma; i++)
		{
			//for (int j = k + 1; j < Na; j++)
			//{
			//	AMatrix[i][j] = AMatrix[i][j] - AMatrix[i][k] * AMatrix[k][j];
			//}
			//AMatrix[i][k] = 0;
			if (AMatrix[i][k] != 0)
			{
				for (int j = k + 1; j < Na; j++)
				{
					AMatrix[i][j] = AMatrix[i][j] - AMatrix[i][k] * AMatrix[k][j];
				}
			}
		}
		k = k + 1;
	}

	double Adiag = AMatrix[Ma - 1][Ma - 1];
	for (int j = Ma - 1; j <= Na - 1; j++)
	{
		AMatrix[Ma - 1][j] = AMatrix[Ma - 1][j] / Adiag;
	}

	int kk = Ma - 1;
	while (kk > 0)
	{
		for (int i = 0; i <= kk - 1; i++)//or for (int i=0;i<kk-1;i++)?
		{
			temp[i] = AMatrix[i][kk];
			for (int j = kk; j <= Na - 1; j++)
			{
				AMatrix[i][j] = AMatrix[i][j] - temp[i] * AMatrix[kk][j];
			}
		}
		kk = kk - 1;
	}


	for (int i = 0; i < Ma; i++)
	{
		I[i] = AMatrix[i][Na - 1];
	}
}

__device__ void solverSMsimple(double AMatrix[Mb][Nb], double I[Mb])
{
	double temp[Mb] = { 0 };
	int k = 0;

	while (k < Mb - 1)
	{
		temp[k] = AMatrix[k][k];

		for (int j = k; j < Nb; j++)
		{
			AMatrix[k][j] = AMatrix[k][j] / temp[k];
		}
		for (int i = k + 1; i < Mb; i++)
		{
			//for (int j = k + 1; j < Na; j++)
			//{
			//	AMatrix[i][j] = AMatrix[i][j] - AMatrix[i][k] * AMatrix[k][j];
			//}
			//AMatrix[i][k] = 0;
			if (AMatrix[i][k] != 0)
			{
				for (int j = k + 1; j < Nb; j++)
				{
					AMatrix[i][j] = AMatrix[i][j] - AMatrix[i][k] * AMatrix[k][j];
				}
			}
		}
		k = k + 1;
	}

	double Adiag = AMatrix[Mb - 1][Mb - 1];
	for (int j = Mb - 1; j <= Nb - 1; j++)
	{
		AMatrix[Mb - 1][j] = AMatrix[Mb - 1][j] / Adiag;
	}

	int kk = Mb - 1;
	while (kk > 0)
	{
		for (int i = 0; i <= kk - 1; i++)//or for (int i=0;i<kk-1;i++)?
		{
			temp[i] = AMatrix[i][kk];
			for (int j = kk; j <= Nb - 1; j++)
			{
				AMatrix[i][j] = AMatrix[i][j] - temp[i] * AMatrix[kk][j];
			}
		}
		kk = kk - 1;
	}


	for (int i = 0; i < Mb; i++)
	{
		I[i] = AMatrix[i][Nb - 1];
	}
}

__device__ void IGBTFWDsimple(double* Geq, double* Ieq, double V_node1, double V_node2, double Vs, double* Itotal)
{
	double IFWD, GFWD, IFWDeq, Rigbt;
	double k = 1.381e-23;
	double q = 1.602e-19;
	double T = 300;
	double VT = k * T / q;
	double M_FWD = 1.21855;
	double ISAT_FWD = 1.001e-6;


	double VFWD = V_node2 - V_node1;



	if (VFWD > 0.7)
	{
		VFWD = 0.7;
	}



	IFWD = ISAT_FWD * (exp(VFWD / M_FWD / VT) - 1);


	GFWD = ISAT_FWD * exp(VFWD / M_FWD / VT) / M_FWD / VT;




	IFWDeq = IFWD - GFWD * VFWD;




	if (Vs > 10)
	{
		Rigbt = 1e-5;
	}

	else
	{
		Rigbt = 1e5;
	}




	*Geq = 1 / Rigbt + GFWD;

	*Ieq = -IFWDeq;

	*Itotal = (V_node1 - V_node2) * *Geq + *Ieq;

}

__device__ void IGBTdiodesimplorer(double G[8][8], double f[8], double V_nodeC, double V_node1, double V_nodeD, double V_nodeG, double V_nodeS, double V_nodeAUX, double V_nodemidFWD,
	double V_nodeE, double h, double Vs, double VGS_old, double QCE_old, double QCG_old, double Qeb_old,
	double* Itotal, double* VGS, double* QCE, double* QCG, double* Qeb)
{
	double VP, K, KLM, A_FET, M_FET, N_FET, BN, M_BJT, ISAT_BJT, Vsat, Isat, VDS, ID;
	double GmosVGS, GmosVDS, Imoseq, ID_Factor1, ID_Factor2, ID_Factor3, d1VDS, d2VDS, d3VDS, d2VGS, d3VGS;
	double k, q, T, VT, VEB, IB, Gbjt, IBeq, RP_BJT, M_FWD, ISAT_FWD, IFWD, GFWD, IFWDeq, VFVD;
	double C_GE, GCGS, ICGS_eq, Rg;
	VP = 6;
	K = 10;
	KLM = 0;
	A_FET = 0.4;
	M_FET = 1.2;
	N_FET = 1.7;
	BN = 15;
	M_BJT = 2;
	ISAT_BJT = 1E-9;


	*VGS = V_nodeG - V_nodeS;

	//2.2 MOSFET
	Vsat = A_FET * pow((*VGS - VP), M_FET);
	Isat = K / 2 * pow((*VGS - VP), N_FET);
	VDS = V_nodeD - V_nodeS;

	if (VDS < 0)
	{
		VDS = 0;
	}

	if (*VGS <= VP)
	{
		ID = 0;

		GmosVGS = 0;
		GmosVDS = 0;
		Imoseq = 0;
	}
	else if (VDS < Vsat)
	{
		ID_Factor1 = 1 + KLM * VDS;
		ID_Factor2 = 2 - VDS / Vsat;
		ID_Factor3 = VDS / Vsat;
		ID = Isat * ID_Factor1 * ID_Factor2 * ID_Factor3;

		d1VDS = KLM;
		d2VDS = -1 / Vsat;
		d3VDS = 1 / Vsat;
		GmosVDS = Isat * (d1VDS * ID_Factor2 * ID_Factor3 + ID_Factor1 * d2VDS * ID_Factor3 + ID_Factor1 * ID_Factor2 * d3VDS);
		d2VGS = VDS / pow(Vsat, 2) * A_FET * M_FET * pow((*VGS - VP), (M_FET - 1));
		d3VGS = -VDS / pow(Vsat, 2) * A_FET * M_FET * pow((*VGS - VP), (M_FET - 1));
		GmosVGS = K / 2 * N_FET * pow((*VGS - VP), (N_FET - 1)) * ID_Factor1 * ID_Factor2 * ID_Factor3 + Isat * ID_Factor1 * d2VGS * ID_Factor3 + Isat * ID_Factor1 * ID_Factor2 * d3VGS;
		Imoseq = ID - GmosVDS * VDS - GmosVGS * (*VGS);
	}
	else
	{
		ID = Isat * (1 + KLM * VDS);

		GmosVDS = Isat * KLM;
		GmosVGS = K / 2 * N_FET * pow((*VGS - VP), (N_FET - 1)) * (1 + KLM * VDS);
		Imoseq = ID - GmosVDS * VDS - GmosVGS * (*VGS);
	}

	//2.3 BJT
	k = 1.381E-23;
	q = 1.602E-19;
	T = 300;
	VT = k * T / q;
	VEB = V_node1 - V_nodeD;

	if (VEB > 2)
	{
		VEB = 2;
	}

	IB = ISAT_BJT * (exp(VEB / M_BJT / VT) - 1);
	Gbjt = ISAT_BJT * exp(VEB / M_BJT / VT) / M_BJT / VT;
	IBeq = IB - Gbjt * VEB;

	RP_BJT = 125;

	////////////2.4 FWD

	M_FWD = 1.21855;
	ISAT_FWD = 1.001e-6;
	VFVD = V_nodeAUX - V_nodemidFWD;
	///////////////////////////////////important////////////////////////////////////
	if (VFVD > 1)
	{
		VFVD = 1;
	}
	////////////////////////////////////////////////////////////////////////////////
	IFWD = ISAT_FWD * (exp(VFVD / M_FWD / VT) - 1);
	GFWD = ISAT_FWD * exp(VFVD / M_FWD / VT) / M_FWD / VT;
	IFWDeq = IFWD - GFWD * VFVD;

	double RB_FWD = 0.0001;
	double INOM = 3000;
	double R_FWD = RB_FWD / sqrt(1 + IFWD / INOM);
	double G_FWD = 1 / R_FWD;

	/////////////2.5 C_GE
	C_GE = 0.7e-7;
	GCGS = C_GE / h;
	ICGS_eq = -GCGS * VGS_old;

	/////////////2.6 C_CE
	double C0_CE = 2e-8;
	double VDIFF_CE = 0.2;
	double VSHIFT_CE = 0;
	double ALPHA_CE = 0.5;
	double BETA_CE = 2;
	double DELTA_CE = 0.0001;
	

	double VCE, B, VVJNCT_CE, CEa, CEb, AA_CE, BB_CE;

	VCE = V_node1 - V_nodeS;    //QCE_old
	B = VSHIFT_CE * VDIFF_CE;
	VVJNCT_CE = B - VCE;
	CEa = BETA_CE - 1;
	CEb = 1 - DELTA_CE;


	AA_CE = (ALPHA_CE * CEb / CEa) / VDIFF_CE;
	BB_CE = (1 - DELTA_CE) * VDIFF_CE / (1 - ALPHA_CE);

	double CCE, ICE, GCE, ICEeq;

	if (VVJNCT_CE >= 0)
	{
		CCE = C0_CE * (1 + CEa * (1 - exp(-VVJNCT_CE * AA_CE)));
		//Charge Calculation by Integration of CCE(from 0 to VCE)
		*QCE = C0_CE * (VCE + CEa * (VCE - exp(VCE * AA_CE) / AA_CE)) + C0_CE * CEa / AA_CE;
		ICE = (*QCE - QCE_old) / h;
		GCE = CCE / h;
		ICEeq = ICE - GCE * VCE;
	}

	else
	{
		CCE = C0_CE * (DELTA_CE + CEb / pow((1 - VVJNCT_CE / VDIFF_CE), ALPHA_CE));
		//Charge Calculation by Integration of CCE(from 0 to VCE)
		*QCE = C0_CE * (DELTA_CE * VCE + BB_CE * pow((1 + VCE / VDIFF_CE), (1 - ALPHA_CE))) - C0_CE * BB_CE;
		ICE = (*QCE - QCE_old) / h;
		GCE = CCE / h;
		ICEeq = ICE - GCE * VCE;
	}

	///////////////////2.7 C_CG
	double C0_CG = 5e-8;
	double VDIFF_CG_OFF = 0.1;
	double VDIFF_CG_ON = 0.1;
	double VSHIFT_CG_OFF = 0;
	double VSHIFT_CG_ON = 0;
	double ALPHA_CG_OFF = 0.5;
	double ALPHA_CG_ON = 0.5;
	double BETA_CG_OFF = 2;
	double BETA_CG_ON = 2;
	double DELTA_CG_OFF = 0.0001;
	double DELTA_CG_ON = 0.0001;

	double VCG, VVJNCT_CG, CGa, CGb, AA_CG, BB_CG, CCG, ICG, GCG, ICGeq;

	VCG = V_node1 - V_nodeG;  //QCG_old
	if (*VGS <= VP)
	{
		B = VSHIFT_CG_OFF * VDIFF_CG_OFF;
		VVJNCT_CG = B - VCG;
		CGa = BETA_CG_OFF - 1;
		CGb = 1 - DELTA_CG_OFF;

		AA_CG = ALPHA_CG_OFF * CGb / CGa / VDIFF_CG_OFF;
		BB_CG = (1 - DELTA_CG_OFF) * VDIFF_CG_OFF / (1 - ALPHA_CG_OFF);

		if (VVJNCT_CG >= 0)
		{
			CCG = C0_CG * (1 + CGa * (1 - exp(-VVJNCT_CG * AA_CG)));
			//Charge Calculation by Integration of CCE(from 0 to VCE)
			*QCG = C0_CG * (VCG + CGa * (VCG - exp(VCG * AA_CG) / AA_CG)) + C0_CG * CGa / AA_CG;
			ICG = (*QCG - QCG_old) / h;
			GCG = CCG / h;
			ICGeq = ICG - GCG * VCG;
		}
		else
		{
			CCG = C0_CG * (DELTA_CG_OFF + CGb / pow((1 - VVJNCT_CG / VDIFF_CG_OFF), ALPHA_CG_OFF));
			//Charge Calculation by Integration of CCE(from 0 to VCE)
			*QCG = C0_CG * (DELTA_CG_OFF * VCG + BB_CG * pow((1 + VCG / VDIFF_CG_OFF), (1 - ALPHA_CG_OFF))) - C0_CG * BB_CG;
			ICG = (*QCG - QCG_old) / h;
			GCG = CCG / h;
			ICGeq = ICG - GCG * VCG;
		}
	}
	else
	{
		B = VSHIFT_CG_ON * VDIFF_CG_ON;
		VVJNCT_CG = B - VCG;
		CGa = BETA_CG_ON - 1;
		CGb = 1 - DELTA_CG_ON;

		AA_CG = ALPHA_CG_ON * CGb / CGa / VDIFF_CG_ON;
		BB_CG = (1 - DELTA_CG_ON) * VDIFF_CG_ON / (1 - ALPHA_CG_ON);

		if (VVJNCT_CG >= 0)
		{
			CCG = C0_CG * (1 + CGa * (1 - exp(-VVJNCT_CG * AA_CG)));
			//Charge Calculation by Integration of CCE(from 0 to VCE)
			*QCG = C0_CG * (VCG + CGa * (VCG - exp(VCG * AA_CG) / AA_CG)) + C0_CG * CGa / AA_CG;
			ICG = (*QCG - QCG_old) / h;
			GCG = CCG / h;
			ICGeq = ICG - GCG * VCG;
		}
		else
		{
			CCG = C0_CG * (DELTA_CG_ON + CGb / pow((1 - VVJNCT_CG / VDIFF_CG_ON), ALPHA_CG_ON));
			//Charge Calculation by Integration of CCE(from 0 to VCE)
			*QCG = C0_CG * (DELTA_CG_ON * VCG + BB_CG * pow((1 + VCG / VDIFF_CG_ON), (1 - ALPHA_CG_ON))) - C0_CG * BB_CG;
			ICG = (*QCG - QCG_old) / h;
			GCG = CCG / h;
			ICGeq = ICG - GCG * VCG;
		}
	}

	//////////////2.8 Ceb(diff)
	double TAU_BE = 5e-8;
	double Ceb, Ieb, GCeb, ICebeq;

	Ceb = TAU_BE * (IB + ISAT_BJT) / M_BJT / VT;
	*Qeb = TAU_BE * ISAT_BJT * (exp(VEB / M_BJT / VT) - 1);
	Ieb = (*Qeb - Qeb_old) / h;
	GCeb = Ceb / h;
	ICebeq = Ieb - GCeb * VEB;


	/////////////2.6 Connector

	Rg = 10.2;
	double LG = 0;
	double LC = 5e-9;
	double RC = 1.15625e-5;
	double LE = 5e-9;
	double RE = 1.15625e-5;
	double RAUX = 1.15625e-5;
	double LAUX = 1.75e-10;


	f[0] = 0;                                             //node C
	f[1] = -(BN + 1) * IBeq - ICEeq - ICGeq - ICebeq;               //node 1
	f[2] = -Imoseq + IBeq + ICebeq;                       //node D  
	f[3] = Vs / Rg - ICGS_eq + ICGeq;                          //node G
	f[4] = Imoseq + BN * IBeq + ICGS_eq + ICEeq;             //node S
	f[5] = -IFWDeq - Vs / Rg;                                //node AUX
	f[6] = IFWDeq;                                        //node mid FWD
	f[7] = 0;                                             //node E

	G[0][0] = 1 / RC;
	G[0][1] = -1 / RC;
	G[0][2] = 0;
	G[0][3] = 0;
	G[0][4] = 0;
	G[0][5] = 0;
	G[0][6] = 0;
	G[0][7] = 0;

	G[1][0] = -1 / RC;
	G[1][1] = (BN + 1) * Gbjt + 1 / RP_BJT + G_FWD + 1 / RC + GCE + GCG + GCeb;
	G[1][2] = -(BN + 1) * Gbjt - 1 / RP_BJT - GCeb;
	G[1][3] = -GCG;
	G[1][4] = -GCE;
	G[1][5] = 0;
	G[1][6] = -G_FWD;
	G[1][7] = 0;

	G[2][0] = 0;
	G[2][1] = -Gbjt - 1 / RP_BJT - GCeb;
	G[2][2] = GmosVDS + Gbjt + 1 / RP_BJT + GCeb;
	G[2][3] = GmosVGS;
	G[2][4] = -GmosVDS - GmosVGS;
	G[2][5] = 0;
	G[2][6] = 0;
	G[2][7] = 0;

	G[3][0] = 0;
	G[3][1] = -GCG;
	G[3][2] = 0;
	G[3][3] = 1 / Rg + GCGS + GCG;
	G[3][4] = -1 / Rg - GCGS;
	G[3][5] = 0;
	G[3][6] = 0;
	G[3][7] = 0;

	G[4][0] = 0;
	G[4][1] = -BN * Gbjt - GCE;
	G[4][2] = -GmosVDS + BN * Gbjt;
	G[4][3] = -GmosVGS - GCGS;
	G[4][4] = GmosVDS + GmosVGS + GCGS + 1 / RAUX + GCE;
	G[4][5] = -1 / RAUX;
	G[4][6] = 0;
	G[4][7] = 0;

	G[5][0] = 0;
	G[5][1] = 0;
	G[5][2] = 0;
	G[5][3] = -1 / Rg;
	G[5][4] = -1 / RAUX;
	G[5][5] = GFWD + 1 / Rg + 1 / RAUX + 1 / RE;
	G[5][6] = -GFWD;
	G[5][7] = -1 / RE;

	G[6][0] = 0;
	G[6][1] = -G_FWD;
	G[6][2] = 0;
	G[6][3] = 0;
	G[6][4] = 0;
	G[6][5] = -GFWD;
	G[6][6] = G_FWD + GFWD;
	G[6][7] = 0;

	G[7][0] = 0;
	G[7][1] = 0;
	G[7][2] = 0;
	G[7][3] = 0;
	G[7][4] = 0;
	G[7][5] = -1 / RE;
	G[7][6] = 0;
	G[7][7] = 1 / RE;



	*Itotal = (V_nodeC - V_node1) / RC;

}

__device__ void submodulesimple(double dt, double U[Mb], double Vg[2], double* Vc, double* vi, double* Vp, double Js,
	double* Ice1, double* Ice2)
{
	double C = 6e-3, ZC = dt / (2 * C);
	double Geq11, Ieq11, Geq22, Ieq22;
	double* Geq1 = &Geq11, * Ieq1 = &Ieq11, * Geq2 = &Geq22, * Ieq2 = &Ieq22;
	double err1 = 1, err1_temp[Mb] = { 1 }, Unew[Mb], U1[2], U2[2];

	for (int i = 0; i < Mb; i++)
	{
		err1_temp[i] = 1;
		Unew[i] = 0;
	}

	int count1 = 0;

	while ((err1 > 1e-6) && (count1 < 50))
	{
		count1 = count1 + 1;

		U1[0] = U[0];
		U1[1] = U[1];
		U2[0] = U[1];
		U2[1] = 0;

		IGBTFWDsimple(Geq1, Ieq1, U1[0], U1[1], Vg[0], Ice1);
		IGBTFWDsimple(Geq2, Ieq2, U2[0], 0, Vg[1], Ice2);

		double Gf[Mb][Nb] = {
			1 / ZC + *Geq1, -*Geq1,-*Ieq1 + 2 * (*vi) / ZC,
			-*Geq1, *Geq1 + *Geq2 ,*Ieq1 - *Ieq2 + Js };
		solverSMsimple(Gf, Unew);

		for (int kk = 0; kk < Mb; kk++)
		{
			err1_temp[kk] = abs((Unew[kk] - U[kk]) / U[kk]);
		}

		err1 = err1_temp[0];
		for (int kk = 0; kk < Mb; kk++)
		{
			if (err1_temp[kk] >= err1)
			{
				err1 = err1_temp[kk];
			}
		}

		for (int kk = 0; kk < Mb; kk++)
		{
			U[kk] = Unew[kk];
		}
	}


	*vi = U1[0] - *vi;
	*Vc = U1[0];
	*Vp = U1[1];



}

__device__ void submodule(double dt, double U[Ma], double Vg[2], double* Vc, double* vi, double* Vp, double Js,
	double VGS_old[2], double QCE_old[2], double QCG_old[2], double Qeb_old[2], double* Ice1, double* Ice2)
{
	//int idx = threadIdx.x+blockIdx.x*blockDim.x;

	double C = 1e-3, ZC = dt / (2 * C);
	double G1[8][8], f1[8], G2[8][8], f2[8], VGS[2], QCE[2], QCG[2], Qeb[2];


	//simulation starts

	double err1 = 1, err1_temp[Ma] = { 1 }, Unew[Ma], U1[8], U2[8];
	for (int i = 0; i < Ma; i++)
	{
		err1_temp[i] = 1;
		Unew[i] = 0;
	}

	int count1 = 0;

	while ((err1 > 1e-6) && (count1 < 30))
	{
		count1 = count1 + 1;


		for (int i = 0; i < 8; i++)
		{
			U1[i] = U[i];
		}
		U2[0] = U[7];
		U2[1] = U[8];
		U2[2] = U[9];
		U2[3] = U[10];
		U2[4] = U[11];
		U2[5] = U[12];
		U2[6] = U[13];
		U2[7] = 0;

		IGBTdiodesimplorer(G1, f1, U1[0], U1[1], U1[2], U1[3], U1[4], U1[5], U1[6], U1[7], dt, Vg[0], VGS_old[0], QCE_old[0], QCG_old[0], Qeb_old[0], Ice1, &VGS[0], &QCE[0], &QCG[0], &Qeb[0]);
		IGBTdiodesimplorer(G2, f2, U2[0], U2[1], U2[2], U2[3], U2[4], U2[5], U2[6], U2[7], dt, Vg[1], VGS_old[1], QCE_old[1], QCG_old[1], Qeb_old[1], Ice2, &VGS[1], &QCE[1], &QCG[1], &Qeb[1]);

		double Gf[Ma][Na] = {
			 1 / ZC + G1[0][0],G1[0][1],G1[0][2],G1[0][3],G1[0][4],G1[0][5],G1[0][6],G1[0][7],0,0,0,0,0,0,f1[0] + 2 * (*vi) / ZC ,
			 G1[1][0],G1[1][1],G1[1][2],G1[1][3],G1[1][4],G1[1][5],G1[1][6],G1[1][7],0,0,0,0,0,0,f1[1] ,
			 G1[2][0],G1[2][1],G1[2][2],G1[2][3],G1[2][4],G1[2][5],G1[2][6],G1[2][7],0,0,0,0,0,0,f1[2] ,
			 G1[3][0],G1[3][1],G1[3][2],G1[3][3],G1[3][4],G1[3][5],G1[3][6],G1[3][7],0,0,0,0,0,0,f1[3] ,
			 G1[4][0],G1[4][1],G1[4][2],G1[4][3],G1[4][4],G1[4][5],G1[4][6],G1[4][7],0,0,0,0,0,0,f1[4] ,
			 G1[5][0],G1[5][1],G1[5][2],G1[5][3],G1[5][4],G1[5][5],G1[5][6],G1[5][7],0,0,0,0,0,0,f1[5] ,
			 G1[6][0],G1[6][1],G1[6][2],G1[6][3],G1[6][4],G1[6][5],G1[6][6],G1[6][7],0,0,0,0,0,0,f1[6] ,
			 G1[7][0],G1[7][1],G1[7][2],G1[7][3],G1[7][4],G1[7][5],G1[7][6],G1[7][7] + 1e-6 + G2[0][0],G2[0][1],G2[0][2],G2[0][3],G2[0][4],G2[0][5],G2[0][6],f1[7] + f2[0] + Js ,
			 0,0,0,0,0,0,0,G2[1][0],G2[1][1],G2[1][2],G2[1][3],G2[1][4],G2[1][5],G2[1][6],f2[1] ,
			 0,0,0,0,0,0,0,G2[2][0],G2[2][1],G2[2][2],G2[2][3],G2[2][4],G2[2][5],G2[2][6],f2[2] ,
			 0,0,0,0,0,0,0,G2[3][0],G2[3][1],G2[3][2],G2[3][3],G2[3][4],G2[3][5],G2[3][6],f2[3] ,
			 0,0,0,0,0,0,0,G2[4][0],G2[4][1],G2[4][2],G2[4][3],G2[4][4],G2[4][5],G2[4][6],f2[4] ,
			 0,0,0,0,0,0,0,G2[5][0],G2[5][1],G2[5][2],G2[5][3],G2[5][4],G2[5][5],G2[5][6],f2[5] ,
			 0,0,0,0,0,0,0,G2[6][0],G2[6][1],G2[6][2],G2[6][3],G2[6][4],G2[6][5],G2[6][6],f2[6] };



		solverSM(Gf, Unew);

		for (int kk = 0; kk < Ma; kk++)
		{
			err1_temp[kk] = abs((Unew[kk] - U[kk]) / U[kk]);
		}

		err1 = err1_temp[0];
		for (int kk = 0; kk < Ma; kk++)
		{
			if (err1_temp[kk] >= err1)
			{
				err1 = err1_temp[kk];
			}
		}

		for (int kk = 0; kk < Ma; kk++)
		{
			U[kk] = Unew[kk];
		}


	}

	////////////////////////////////////////////*Ice1 = vcge_old[0];//2*(*vi)/ZC-U1[0]/ZC;
	////////////////////////////////////////////*Ice2 = vcge[0];//*Ice1+Js;


	for (int i = 0; i < 2; i++)// move to the outer layer??
	{
		VGS_old[i] = VGS[i];
		QCE_old[i] = QCE[i];
		QCG_old[i] = QCG[i];
		Qeb_old[i] = Qeb[i];
	}


	*vi = U1[0] - *vi;
	*Vc = U1[0];
	*Vp = U1[7];

	//__syncthreads();

}


__global__ void PQcontrol(double* theta, double* Uref, double* sum1, double* sum2, double* Idaxis, double* Iqaxis, double* sum3, double* sum4, double* sum5, double* Vdc)
{
	int idx = threadIdx.x + blockIdx.x * blockDim.x;
	double dt = 1e-8;


	double Vdc_ref = 60e3;

	double Vsd = 31e3 / sqrtf(3) * sqrtf(2), Vsq = 0;

	Vsd = Vsd / (31e3 / sqrtf(3) * sqrtf(2));
	double Pref = 1, Qref = 0;
	//double evdc, Idref, ed, ud, md, Kvs = 60000/1e3, Kvp = 10000/1e3, Ki = 500000 / 1e3, Kp = 10000 / 1e3, Vsq = 0, L = 10e-3,
	//	eq, uq, mq, Iq_ref = 0, alpha, beta;
	double evdc, ep, Idref, ed, ud, md, Kvs = 100, Ki = 6, Kp = 0.6, Larm_pu = 0.15, Rarm_pu = Larm_pu / 100,
		eq, erq, uq, mq, Iq_ref = 0, alpha, beta;
	double Kvp = 4;
	//double Kvp = 40;
	double Lff = Larm_pu * 0.4 + Larm_pu / 2, Rff = Rarm_pu / 2;
	double Kp_Preg = 0.5 / 3, Ki_Preg = 1.0;
	double Kp_Qreg = 0.5 / 3, Ki_Qreg = 1.0;
	double Pmeas = Vsd * Idaxis[idx];
	double Qmeas = -Vsd * Iqaxis[idx];
	//dq control
	theta[idx] = theta[idx] + 2 * PI * freq * dt;
	__syncthreads();

	double Vnom_dc = 60e3;

	double Vnom_sec = 31e3;

	if (idx == 0)//0, rectifier
	{
		ep = Pmeas - Pref;
		sum4[idx] = sum4[idx] + Ki_Preg * dt * ep;
		Vdc_ref = ep * Kp_Preg + sum4[idx];
		if (Vdc_ref > 1.2)
		{
			Vdc_ref = 1.2;
		}
		if (Vdc_ref < 0.8)
		{
			Vdc_ref = 0.8;
		}
		Vdc_ref = Vdc_ref * Vnom_dc;
		evdc = Vdc[idx] - Vdc_ref;
		evdc = evdc / Vnom_dc;
		sum1[idx] = sum1[idx] + Kvs * dt * evdc;
		Idref = evdc * Kvp + sum1[idx];
	}
	else//1, inverter
	{
		evdc = Vdc[idx] - Vdc_ref;
		evdc = evdc / Vnom_dc;
		sum1[idx] = sum1[idx] + Kvs * dt * evdc;
		Idref = evdc * Kvp + sum1[idx];
	}
	if (Idref > 2)
	{
		Idref = 2;
	}
	if (Idref < -2)
	{
		Idref = -2;
	}
	__syncthreads();




	erq = Qref - Qmeas;
	sum5[idx] = sum5[idx] + Ki_Qreg * dt * erq;
	Iq_ref = erq * Kp_Qreg + sum5[idx];
	if (Iq_ref > 0.25)
	{
		Iq_ref = 0.25;
	}
	if (Iq_ref < -0.25)
	{
		Iq_ref = -0.25;
	}
	Iq_ref = -Iq_ref - Qref / Vsd;
	if (Iq_ref > 0.25)
	{
		Iq_ref = 0.25;
	}
	if (Iq_ref < -0.25)
	{
		Iq_ref = -0.25;
	}

	ed = Idref - Idaxis[idx];
	sum2[idx] = sum2[idx] + Ki * dt * ed;
	ud = sum2[idx] + Kp * ed;
	if (ud > 2)
	{
		ud = 2;
	}
	if (ud < -2)
	{
		ud = -2;
	}
	double Pnom = 18e6;

	md = ud + Vsd + Idref * Rff - Lff * Iq_ref;

	eq = Iq_ref - Iqaxis[idx];
	sum3[idx] = sum3[idx] + Ki * dt * eq;
	uq = sum3[idx] + Kp * eq;
	if (uq > 2)
	{
		uq = 2;
	}
	if (uq < -2)
	{
		uq = -2;
	}
	mq = uq + Vsq + Iq_ref * Rff + Lff * Idref;

	double Vdc_val = Vdc[idx];
	if (Vdc_val > Vnom_dc * 1.25)
	{
		Vdc_val = Vnom_dc * 1.25;
	}
	if (Vdc_val < Vnom_dc * 0.75)
	{
		Vdc_val = Vnom_dc * 0.75;
	}

	md = md * Vnom_dc / Vdc_val;
	mq = mq * Vnom_dc / Vdc_val;

	md = md / (Vnom_dc / 2 * sqrtf(3) / sqrtf(2) / Vnom_sec);
	mq = mq / (Vnom_dc / 2 * sqrtf(3) / sqrtf(2) / Vnom_sec);
	if (mq > 2)
	{
		mq = 2;
	}
	if (mq < -2)
	{
		mq = -2;
	}
	if (md > 2)
	{
		md = 2;
	}
	if (md < -2)
	{
		md = -2;
	}
	alpha = md * cosf(theta[idx]) - mq * sinf(theta[idx]);
	beta = md * sinf(theta[idx]) + mq * cosf(theta[idx]);
	Uref[3 * idx] = alpha * (Vnom_sec / sqrtf(3) * sqrtf(2));
	Uref[1 + 3 * idx] = (-0.5 * alpha + sqrt(3.0) * beta / 2) * (Vnom_sec / sqrtf(3) * sqrtf(2));
	Uref[2 + 3 * idx] = (-0.5 * alpha - sqrt(3.0) * beta / 2) * (Vnom_sec / sqrtf(3) * sqrtf(2));


	__syncthreads();
}

__global__ void PQcontrol2(double* theta, double* Uref, double* sum1, double* sum2, double* Idaxis, double* Iqaxis, double* sum3, double* sum4, double* sum5, double* Vdc)

{
	int idx = threadIdx.x + blockIdx.x * blockDim.x;
	double dt = 1e-8;


	double Vdc_ref = 60e3;

	double Vsd = 31e3 / sqrtf(3) * sqrtf(2), Vsq = 0;

	Vsd = Vsd / (31e3 / sqrtf(3) * sqrtf(2));
	double Pref = 1, Qref = 0;
	//double evdc, Idref, ed, ud, md, Kvs = 60000/1e3, Kvp = 10000/1e3, Ki = 500000 / 1e3, Kp = 10000 / 1e3, Vsq = 0, L = 10e-3,
	//	eq, uq, mq, Iq_ref = 0, alpha, beta;
	double evdc, ep, Idref, ed, ud, md, Kvs = 100, Ki = 6, Kp = 0.6, Larm_pu = 0.15, Rarm_pu = Larm_pu / 100,
		eq, erq, uq, mq, Iq_ref = 0, alpha, beta;
	double Kvp = 4;
	//double Kvp = 40;
	double Lff = Larm_pu * 0.4 + Larm_pu / 2, Rff = Rarm_pu / 2;
	double Kp_Preg = 0.5 / 3, Ki_Preg = 1.0;
	double Kp_Qreg = 0.5 / 3, Ki_Qreg = 1.0;
	double Pmeas = Vsd * Idaxis[idx];
	double Qmeas = -Vsd * Iqaxis[idx];
	//dq control
	theta[idx] = theta[idx] + 2 * PI * freq * dt;
	__syncthreads();

	double Vnom_dc = 60e3;

	double Vnom_sec = 31e3;

	if (idx == 0)//0, rectifier
	{
		ep = Pmeas - Pref;
		sum4[idx] = sum4[idx] + Ki_Preg * dt * ep;
		Vdc_ref = ep * Kp_Preg + sum4[idx];
		if (Vdc_ref > 1.2)
		{
			Vdc_ref = 1.2;
		}
		if (Vdc_ref < 0.8)
		{
			Vdc_ref = 0.8;
		}
		Vdc_ref = Vdc_ref * Vnom_dc;
		evdc = Vdc[idx] - Vdc_ref;
		evdc = evdc / Vnom_dc;
		sum1[idx] = sum1[idx] + Kvs * dt * evdc;
		Idref = evdc * Kvp + sum1[idx];
	}
	else//1, inverter
	{
		evdc = Vdc[idx] - Vdc_ref;
		evdc = evdc / Vnom_dc;
		sum1[idx] = sum1[idx] + Kvs * dt * evdc;
		Idref = evdc * Kvp + sum1[idx];
	}
	if (Idref > 2)
	{
		Idref = 2;
	}
	if (Idref < -2)
	{
		Idref = -2;
	}
	__syncthreads();




	erq = Qref - Qmeas;
	sum5[idx] = sum5[idx] + Ki_Qreg * dt * erq;
	Iq_ref = erq * Kp_Qreg + sum5[idx];
	if (Iq_ref > 0.25)
	{
		Iq_ref = 0.25;
	}
	if (Iq_ref < -0.25)
	{
		Iq_ref = -0.25;
	}
	Iq_ref = -Iq_ref - Qref / Vsd;
	if (Iq_ref > 0.25)
	{
		Iq_ref = 0.25;
	}
	if (Iq_ref < -0.25)
	{
		Iq_ref = -0.25;
	}

	ed = Idref - Idaxis[idx];
	sum2[idx] = sum2[idx] + Ki * dt * ed;
	ud = sum2[idx] + Kp * ed;
	if (ud > 2)
	{
		ud = 2;
	}
	if (ud < -2)
	{
		ud = -2;
	}
	double Pnom = 18e6;

	md = ud + Vsd + Idref * Rff - Lff * Iq_ref;

	eq = Iq_ref - Iqaxis[idx];
	sum3[idx] = sum3[idx] + Ki * dt * eq;
	uq = sum3[idx] + Kp * eq;
	if (uq > 2)
	{
		uq = 2;
	}
	if (uq < -2)
	{
		uq = -2;
	}
	mq = uq + Vsq + Iq_ref * Rff + Lff * Idref;

	double Vdc_val = Vdc[idx];
	if (Vdc_val > Vnom_dc * 1.25)
	{
		Vdc_val = Vnom_dc * 1.25;
	}
	if (Vdc_val < Vnom_dc * 0.75)
	{
		Vdc_val = Vnom_dc * 0.75;
	}

	md = md * Vnom_dc / Vdc_val;
	mq = mq * Vnom_dc / Vdc_val;

	md = md / (Vnom_dc / 2 * sqrtf(3) / sqrtf(2) / Vnom_sec);
	mq = mq / (Vnom_dc / 2 * sqrtf(3) / sqrtf(2) / Vnom_sec);
	if (mq > 2)
	{
		mq = 2;
	}
	if (mq < -2)
	{
		mq = -2;
	}
	if (md > 2)
	{
		md = 2;
	}
	if (md < -2)
	{
		md = -2;
	}
	alpha = md * cosf(theta[idx]) - mq * sinf(theta[idx]);
	beta = md * sinf(theta[idx]) + mq * cosf(theta[idx]);
	Uref[3 * idx] = alpha * (Vnom_sec / sqrtf(3) * sqrtf(2));
	Uref[1 + 3 * idx] = (-0.5 * alpha + sqrt(3.0) * beta / 2) * (Vnom_sec / sqrtf(3) * sqrtf(2));
	Uref[2 + 3 * idx] = (-0.5 * alpha - sqrt(3.0) * beta / 2) * (Vnom_sec / sqrtf(3) * sqrtf(2));


	__syncthreads();
}

__global__ void PSC(double* Vc, double* Vzu2, double* Vau2, double* Vau, double* Vdc, double* Js1, double* Js2)
{
	int idx = threadIdx.x + blockIdx.x * blockDim.x;//leg No. 6: 0,1,2,3,4,5

	double dt = 1e-5;
	double Vcref, K1 = 0.5 / 10, K2 = 1 / 10, K3 = 0.5 / 10, K4 = 1 / 10;

	double Vzu1, Vzu, Vzuref, Vau1;//one phase
	double Vcsum;
	double Iu, Id, Vcave;

	//-----------------------------------------leg No. 6
	Vcsum = 0;
	for (int k = idx * NoM; k < idx * NoM + NoM; k++)
	{
		Vcsum = Vcsum + Vc[k];
	}


	Vcave = Vcsum / NoM;

	int jmmc = idx / 3;//0,0,0,1,1,1

	Vcref = 1.0 * Vdc[jmmc] / NoC;//Vcref=1.0*200.0/NoC;//

	Iu = Js1[idx];//0,1,2,3,4,5
	Id = Js2[idx];//0,1,2,3,4,5
	Vzu1 = (Vcref - Vcave) * K1;
	Vzu2[idx] = Vzu2[idx] + (Vcref - Vcave) * dt * K2;
	Vzu = Vzu1 + Vzu2[idx];
	Vzuref = (Iu + Id) / 2;
	Vau1 = (Vzuref - Vzu) * K3;
	Vau2[idx] = Vau2[idx] + (Vzuref - Vzu) * dt * K4;
	Vau[idx] = Vau1 + Vau2[idx];

	__syncthreads();

}

__global__ void NLM(int t, double* Vc, double* Uref, double* Js1, double* Js2, double* Vg, double* Deadtime_count, double* NNU, double* NND)
{
	int idx = threadIdx.x + blockIdx.x * blockDim.x;//leg No. 6: 0,1,2,3,4,5

	double dt = 1e-5;
	//Tctr = 100e-6;
	if (t % 1000 == 0)
	{
		int NU, ND;

		//-----------------------------------------leg No. 6
		NU = NoC / 2 - round(Uref[idx] / (200e3 / NoC));
		ND = NoC / 2 + round(Uref[idx] / (200e3 / NoC));
		if (NU < 0)
		{
			NU = 0;
		}
		if (NU > NoC)
		{
			NU = NoC;
		}
		if (ND < 0)
		{
			ND = 0;
		}
		if (ND > NoC)
		{
			ND = NoC;
		}
		*NNU = NU;
		*NND = ND;
		//NU = roundf(Urefp[idx] / (200 / NoC));
		//ND = roundf(Urefn[idx] / (200 / NoC));

		//////////////////////Number for sub modules from up to down//////////////////////
		int idxU1[NoC], idxU2[NoC];

		for (int i = 0; i < NoC; i++)
		{
			idxU1[i] = idx * NoM + i;
			idxU2[i] = idx * NoM + i + NoC;
		}


		/////////////////////////////////////////////////////////////////////////////////


		///////////////////////extract capacitance voltage for all sub modules///////////

		double Vcu[NoC], Vcd[NoC];
		for (int i = 0; i < NoC; i++)
		{
			Vcu[i] = Vc[idx * NoM + i];
			Vcd[i] = Vc[idx * NoM + i + NoC];
		}

		for (int i = 0; i < NoC - 1; i++)
		{
			for (int j = 0; j < NoC - i - 1; j++)
			{
				if (Vcu[j] > Vcu[j + 1])
				{
					double temp1 = Vcu[j]; Vcu[j] = Vcu[j + 1]; Vcu[j + 1] = temp1;
					double temp2 = idxU1[j]; idxU1[j] = idxU1[j + 1]; idxU1[j + 1] = temp2;
				}
				if (Vcd[j] > Vcd[j + 1])
				{
					double temp1 = Vcd[j]; Vcd[j] = Vcd[j + 1]; Vcd[j + 1] = temp1;
					double temp2 = idxU2[j]; idxU2[j] = idxU2[j + 1]; idxU2[j + 1] = temp2;
				}
			}
		}


		int k1, k2;
		if (Js1[idx] < 0)
		{
			for (int i = 0; i < NU; i++)
			{
				k1 = idxU1[NoC - 1 - i];
				k2 = 2 * k1;
				Vg[k2] = 15;
				//Vg[k2] = -10;
				//Deadtime_count[k2] = count;
				Vg[k2 + 1] = -10;
			}
			for (int i = NU; i < NoC; i++)

			{
				k1 = idxU1[NoC - 1 - i];
				k2 = 2 * k1;
				Vg[k2] = -10;
				Vg[k2 + 1] = 15;
				//Vg[k2 + 1] = -10;
				//Deadtime_count[k2 + 1] = count;
			}
		}

		if (Js1[idx] >= 0)
		{
			for (int i = 0; i < NU; i++)
			{
				k1 = idxU1[i];
				k2 = 2 * k1;
				Vg[k2] = 15;
				//Vg[k2] = -10;
				//Deadtime_count[k2] = count;
				Vg[k2 + 1] = -10;
			}
			for (int i = NU; i < NoC; i++)

			{
				k1 = idxU1[i];
				k2 = 2 * k1;
				Vg[k2] = -10;
				Vg[k2 + 1] = 15;
				//Vg[k2 + 1] = -10;
				//Deadtime_count[k2 + 1] = count;
			}
		}



		if (Js2[idx] < 0)
		{
			for (int i = 0; i < ND; i++)

			{
				k1 = idxU2[NoC - 1 - i];
				k2 = 2 * k1;
				Vg[k2] = 15;
				//Vg[k2] = -10;
				//Deadtime_count[k2] = count;
				Vg[k2 + 1] = -10;
			}
			for (int i = ND; i < NoC; i++)

			{
				k1 = idxU2[NoC - 1 - i];
				k2 = 2 * k1;
				Vg[k2] = -10;
				Vg[k2 + 1] = 15;
				//Vg[k2 + 1] = -10;
				//Deadtime_count[k2 + 1] = count;
			}
		}

		if (Js2[idx] >= 0)
		{
			for (int i = 0; i < ND; i++)

			{
				k1 = idxU2[i];
				k2 = 2 * k1;
				Vg[k2] = 15;
				//Vg[k2] = -10;
				//Deadtime_count[k2] = count;
				Vg[k2 + 1] = -10;
			}
			for (int i = ND; i < NoC; i++)

			{
				k1 = idxU2[i];
				k2 = 2 * k1;
				Vg[k2] = -10;
				Vg[k2 + 1] = 15;
				//Vg[k2 + 1] = -10;
				//Deadtime_count[k2 + 1] = count;
			}
		}
		//////////////////////////////////////////////////////////////////////////////////////////////////
		/*
		for (int i = 0; i < NU; i++)
		{
			Vg[2 * (idx * NoM + i)] = 15;
			//Vg[2 * (idx * NoM + i)] = -10;
			//Deadtime_count[2 * (idx * NoM + i)] = count;
			Vg[2 * (idx * NoM + i) + 1] = -10;
		}
		for (int i = NU; i < NoC; i++)
		{
			Vg[2 * (idx * NoM + i)] = -10;
			Vg[2 * (idx * NoM + i) + 1] = 15;
			//Vg[2 * (idx * NoM + i) + 1] = -10;
			//Deadtime_count[2 * (idx * NoM + i) + 1] = count;
		}

		for (int i = NoC; i < NoC + ND; i++)
		{
			Vg[2 * (idx * NoM + i)] = 15;
			//Vg[2 * (idx * NoM + i)] = -10;
			//Deadtime_count[2 * (idx * NoM + i)] = count;
			Vg[2 * (idx * NoM + i) + 1] = -10;
		}
		for (int i = NoC + ND; i < NoM; i++)
		{
			Vg[2 * (idx * NoM + i)] = -10;
			Vg[2 * (idx * NoM + i) + 1] = 15;
			//Vg[2 * (idx * NoM + i) + 1] = -10;
			//Deadtime_count[2 * (idx * NoM + i) + 1] = count;
		}
		*/






	}


	//for (int i = 2 * idx * NoM; i < 2 * idx * NoM + 2 * NoM; i++)

	//{
	//	if (Deadtime_count[i] == 0)
	//	{
	//		Vg[i] = 15;
	//	}
	//	Deadtime_count[i] = Deadtime_count[i] - 1;
	//}


	__syncthreads();

}

//void NLM2(int t, double *Vc, double *Uref, double *Js1, double *Js2, double *Vg, double *Deadtime_count, double *NNU, double *NND)
__global__ void NLM2(int t, double* Vc, double* Uref, double* Js1, double* Js2, double* Vg, double* Deadtime_count, double* NNU, double* NND)
{
	int idx = threadIdx.x + blockIdx.x * blockDim.x;//leg No. 6: 0,1,2,3,4,5
	//for (int idx = 0; idx < 6; idx++)
	//{
	double dt = 1e-8;
	//Tctr = 100e-6;

	int NU, ND;

	//-----------------------------------------leg No. 6
	NU = NoC1 / 2 - round(Uref[idx] / (60e3 / NoC1));
	ND = NoC1 / 2 + round(Uref[idx] / (60e3 / NoC1));
	if (NU < 0)
	{
		NU = 0;
	}
	if (NU > NoC1)
	{
		NU = NoC1;
	}
	if (ND < 0)
	{
		ND = 0;
	}
	if (ND > NoC1)
	{
		ND = NoC1;
	}

	if (t % 10000 == 0)
	{

		NNU[idx] = NU;
		NND[idx] = ND;
		//NU = roundf(Urefp[idx] / (200 / NoC));
		//ND = roundf(Urefn[idx] / (200 / NoC));

		//////////////////////Number for sub modules from up to down//////////////////////
		int idxU1[NoC], idxU2[NoC];

		for (int i = 0; i < NoC; i++)
		{
			idxU1[i] = i;
			idxU2[i] = i;
		}


		/////////////////////////////////////////////////////////////////////////////////


		///////////////////////extract capacitance voltage for all sub modules///////////

		double Vcu[NoC], Vcd[NoC];
		for (int i = 0; i < NoC; i++)
		{
			Vcu[i] = Vc[idx * NoM + i];
			Vcd[i] = Vc[idx * NoM + i + NoC];
		}

		for (int i = 0; i < NoC - 1; i++)
		{
			for (int j = 0; j < NoC - i - 1; j++)
			{
				if (Vcu[j] > Vcu[j + 1])
				{
					double temp1 = Vcu[j]; Vcu[j] = Vcu[j + 1]; Vcu[j + 1] = temp1;
					double temp2 = idxU1[j]; idxU1[j] = idxU1[j + 1]; idxU1[j + 1] = temp2;
				}
				if (Vcd[j] > Vcd[j + 1])
				{
					double temp1 = Vcd[j]; Vcd[j] = Vcd[j + 1]; Vcd[j + 1] = temp1;
					double temp2 = idxU2[j]; idxU2[j] = idxU2[j + 1]; idxU2[j + 1] = temp2;
				}
			}
		}
		__syncthreads();

		int k1, k2;
		if (Js1[idx] >= 0)
		{
			for (int i = 0; i < NU; i++)
			{
				k1 = idx * NoM + idxU1[i];
				k2 = 2 * k1;
				Vg[k2] = 15;
				//Vg[k2] = -10;
				//Deadtime_count[k2] = count;
				Vg[k2 + 1] = -10;
			}
			for (int i = NU; i < NoC; i++)

			{
				k1 = idx * NoM + idxU1[i];
				k2 = 2 * k1;
				Vg[k2] = -10;
				Vg[k2 + 1] = 15;
				//Vg[k2 + 1] = -10;
				//Deadtime_count[k2 + 1] = count;
			}
		}
		else
		{

			for (int i = 0; i < NU; i++)
			{
				k1 = idx * NoM + idxU1[NoC - 1 - i];
				k2 = 2 * k1;
				Vg[k2] = 15;
				//Vg[k2] = -10;
				//Deadtime_count[k2] = count;
				Vg[k2 + 1] = -10;
			}
			for (int i = NU; i < NoC; i++)

			{
				k1 = idx * NoM + idxU1[NoC - 1 - i];
				k2 = 2 * k1;
				Vg[k2] = -10;
				Vg[k2 + 1] = 15;
				//Vg[k2 + 1] = -10;
				//Deadtime_count[k2 + 1] = count;
			}
		}

		__syncthreads();

		if (Js2[idx] >= 0)
		{
			for (int i = 0; i < ND; i++)

			{
				k1 = idx * NoM + idxU2[i] + NoC;
				k2 = 2 * k1;
				Vg[k2] = 15;
				//Vg[k2] = -10;
				//Deadtime_count[k2] = count;
				Vg[k2 + 1] = -10;
			}
			for (int i = ND; i < NoC; i++)

			{
				k1 = idx * NoM + idxU2[i] + NoC;
				k2 = 2 * k1;
				Vg[k2] = -10;
				Vg[k2 + 1] = 15;
				//Vg[k2 + 1] = -10;
				//Deadtime_count[k2 + 1] = count;
			}

		}
		else
		{
			for (int i = 0; i < ND; i++)

			{
				k1 = idx * NoM + idxU2[NoC - 1 - i] + NoC;
				k2 = 2 * k1;
				Vg[k2] = 15;
				//Vg[k2] = -10;
				//Deadtime_count[k2] = count;
				Vg[k2 + 1] = -10;
			}
			for (int i = ND; i < NoC; i++)

			{
				k1 = idx * NoM + idxU2[NoC - 1 - i] + NoC;
				k2 = 2 * k1;
				Vg[k2] = -10;
				Vg[k2 + 1] = 15;
				//Vg[k2 + 1] = -10;
				//Deadtime_count[k2 + 1] = count;
			}
		}
		//////////////////////////////////////////////////////////////////////////////////////////////////
		/*
		for (int i = 0; i < NU; i++)
		{
			Vg[2 * (idx * NoM + i)] = 15;
			//Vg[2 * (idx * NoM + i)] = -10;
			//Deadtime_count[2 * (idx * NoM + i)] = count;
			Vg[2 * (idx * NoM + i) + 1] = -10;
		}
		for (int i = NU; i < NoC; i++)
		{
			Vg[2 * (idx * NoM + i)] = -10;
			Vg[2 * (idx * NoM + i) + 1] = 15;
			//Vg[2 * (idx * NoM + i) + 1] = -10;
			//Deadtime_count[2 * (idx * NoM + i) + 1] = count;
		}

		for (int i = NoC; i < NoC + ND; i++)
		{
			Vg[2 * (idx * NoM + i)] = 15;
			//Vg[2 * (idx * NoM + i)] = -10;
			//Deadtime_count[2 * (idx * NoM + i)] = count;
			Vg[2 * (idx * NoM + i) + 1] = -10;
		}
		for (int i = NoC + ND; i < NoM; i++)
		{
			Vg[2 * (idx * NoM + i)] = -10;
			Vg[2 * (idx * NoM + i) + 1] = 15;
			//Vg[2 * (idx * NoM + i) + 1] = -10;
			//Deadtime_count[2 * (idx * NoM + i) + 1] = count;
		}
		*/

	}
	//}


	__syncthreads();

}

__global__ void NLM22(int t, double* Vc, double* Uref, double* Js1, double* Js2, double* Vg, double* Deadtime_count, double* NNU, double* NND)
{
	int idx = threadIdx.x + blockIdx.x * blockDim.x;//leg No. 6: 0,1,2,3,4,5
	//for (int idx = 0; idx < 6; idx++)
	//{
	double dt = 1e-8;
	//Tctr = 100e-6;

	int NU, ND;

	//-----------------------------------------leg No. 6
	NU = NoC1 / 2 - round(Uref[idx] / (60e3 / NoC1));
	ND = NoC1 / 2 + round(Uref[idx] / (60e3 / NoC1));
	if (NU < 0)
	{
		NU = 0;
	}
	if (NU > NoC1)
	{
		NU = NoC1;
	}
	if (ND < 0)
	{
		ND = 0;
	}
	if (ND > NoC1)
	{
		ND = NoC1;
	}

	if (t % 10000 == 0)
	{

		NNU[idx] = NU;
		NND[idx] = ND;
		//NU = roundf(Urefp[idx] / (200 / NoC));
		//ND = roundf(Urefn[idx] / (200 / NoC));

		//////////////////////Number for sub modules from up to down//////////////////////
		int idxU1[NoC], idxU2[NoC];

		for (int i = 0; i < NoC; i++)
		{
			idxU1[i] = i;
			idxU2[i] = i;
		}


		/////////////////////////////////////////////////////////////////////////////////


		///////////////////////extract capacitance voltage for all sub modules///////////

		double Vcu[NoC], Vcd[NoC];
		for (int i = 0; i < NoC; i++)
		{
			Vcu[i] = Vc[idx * NoM + i];
			Vcd[i] = Vc[idx * NoM + i + NoC];
		}

		for (int i = 0; i < NoC - 1; i++)
		{
			for (int j = 0; j < NoC - i - 1; j++)
			{
				if (Vcu[j] > Vcu[j + 1])
				{
					double temp1 = Vcu[j]; Vcu[j] = Vcu[j + 1]; Vcu[j + 1] = temp1;
					double temp2 = idxU1[j]; idxU1[j] = idxU1[j + 1]; idxU1[j + 1] = temp2;
				}
				if (Vcd[j] > Vcd[j + 1])
				{
					double temp1 = Vcd[j]; Vcd[j] = Vcd[j + 1]; Vcd[j + 1] = temp1;
					double temp2 = idxU2[j]; idxU2[j] = idxU2[j + 1]; idxU2[j + 1] = temp2;
				}
			}
		}
		__syncthreads();

		int k1, k2;
		if (Js1[idx] >= 0)
		{
			for (int i = 0; i < NU; i++)
			{
				k1 = idx * NoM + idxU1[i];
				k2 = 2 * k1;
				Vg[k2] = 15;
				//Vg[k2] = -10;
				//Deadtime_count[k2] = count;
				Vg[k2 + 1] = -10;
			}
			for (int i = NU; i < NoC; i++)

			{
				k1 = idx * NoM + idxU1[i];
				k2 = 2 * k1;
				Vg[k2] = -10;
				Vg[k2 + 1] = 15;
				//Vg[k2 + 1] = -10;
				//Deadtime_count[k2 + 1] = count;
			}
		}
		else
		{

			for (int i = 0; i < NU; i++)
			{
				k1 = idx * NoM + idxU1[NoC - 1 - i];
				k2 = 2 * k1;
				Vg[k2] = 15;
				//Vg[k2] = -10;
				//Deadtime_count[k2] = count;
				Vg[k2 + 1] = -10;
			}
			for (int i = NU; i < NoC; i++)

			{
				k1 = idx * NoM + idxU1[NoC - 1 - i];
				k2 = 2 * k1;
				Vg[k2] = -10;
				Vg[k2 + 1] = 15;
				//Vg[k2 + 1] = -10;
				//Deadtime_count[k2 + 1] = count;
			}
		}

		__syncthreads();

		if (Js2[idx] >= 0)
		{
			for (int i = 0; i < ND; i++)

			{
				k1 = idx * NoM + idxU2[i] + NoC;
				k2 = 2 * k1;
				Vg[k2] = 15;
				//Vg[k2] = -10;
				//Deadtime_count[k2] = count;
				Vg[k2 + 1] = -10;
			}
			for (int i = ND; i < NoC; i++)

			{
				k1 = idx * NoM + idxU2[i] + NoC;
				k2 = 2 * k1;
				Vg[k2] = -10;
				Vg[k2 + 1] = 15;
				//Vg[k2 + 1] = -10;
				//Deadtime_count[k2 + 1] = count;
			}

		}
		else
		{
			for (int i = 0; i < ND; i++)

			{
				k1 = idx * NoM + idxU2[NoC - 1 - i] + NoC;
				k2 = 2 * k1;
				Vg[k2] = 15;
				//Vg[k2] = -10;
				//Deadtime_count[k2] = count;
				Vg[k2 + 1] = -10;
			}
			for (int i = ND; i < NoC; i++)

			{
				k1 = idx * NoM + idxU2[NoC - 1 - i] + NoC;
				k2 = 2 * k1;
				Vg[k2] = -10;
				Vg[k2 + 1] = 15;
				//Vg[k2 + 1] = -10;
				//Deadtime_count[k2 + 1] = count;
			}
		}
		//////////////////////////////////////////////////////////////////////////////////////////////////
		/*
		for (int i = 0; i < NU; i++)
		{
			Vg[2 * (idx * NoM + i)] = 15;
			//Vg[2 * (idx * NoM + i)] = -10;
			//Deadtime_count[2 * (idx * NoM + i)] = count;
			Vg[2 * (idx * NoM + i) + 1] = -10;
		}
		for (int i = NU; i < NoC; i++)
		{
			Vg[2 * (idx * NoM + i)] = -10;
			Vg[2 * (idx * NoM + i) + 1] = 15;
			//Vg[2 * (idx * NoM + i) + 1] = -10;
			//Deadtime_count[2 * (idx * NoM + i) + 1] = count;
		}

		for (int i = NoC; i < NoC + ND; i++)
		{
			Vg[2 * (idx * NoM + i)] = 15;
			//Vg[2 * (idx * NoM + i)] = -10;
			//Deadtime_count[2 * (idx * NoM + i)] = count;
			Vg[2 * (idx * NoM + i) + 1] = -10;
		}
		for (int i = NoC + ND; i < NoM; i++)
		{
			Vg[2 * (idx * NoM + i)] = -10;
			Vg[2 * (idx * NoM + i) + 1] = 15;
			//Vg[2 * (idx * NoM + i) + 1] = -10;
			//Deadtime_count[2 * (idx * NoM + i) + 1] = count;
		}
		*/

	}
	//}


	__syncthreads();

}

__global__ void SM(int t, double* Uref, double* Js1, double* Js2, double* vi, double* Vp, double* Vc, double* Vdc,
	double* U, double* VGS_old, double* QCE_old, double* QCG_old, double* Qeb_old, double* Ice1, double* Ice2, double* out1, double* Vgg)
{//48D
	int j = threadIdx.x + blockIdx.x * blockDim.x;


	double dt = 1e-8;


	double Vg[2];
	double Js;

	Vg[0] = Vgg[2 * j];
	Vg[1] = Vgg[2 * j + 1];

	int jleg = j / NoM;

	int jvc = j % NoM;//0,1,2,3,4,5,6,7;0,1,2,3,4,5,6,7;...



	if (jvc < NoC)
	{
		Js = Js1[jleg];
	}
	else
	{
		Js = Js2[jleg];
	}


	if (j == 3)
	{
		*out1 = Vg[1];
	}

	double U_t[Ma], VGS_old_t[2], QCE_old_t[2], QCG_old_t[2], Qeb_old_t[2];
	for (int i = 0; i < Ma; i++)
	{
		U_t[i] = U[j * Ma + i];
	}
	for (int i = 0; i < 2; i++)
	{
		VGS_old_t[i] = VGS_old[j * 2 + i];
		QCE_old_t[i] = QCE_old[j * 2 + i];
		QCG_old_t[i] = QCG_old[j * 2 + i];
		Qeb_old_t[i] = Qeb_old[j * 2 + i];
	}

	submodule(dt, U_t, Vg, &Vc[j], &vi[j], &Vp[j], Js, VGS_old_t, QCE_old_t, QCG_old_t, Qeb_old_t, &Ice1[j], &Ice2[j]);

	for (int i = 0; i < Ma; i++)
	{
		U[j * Ma + i] = U_t[i];
	}
	for (int i = 0; i < 2; i++)
	{
		VGS_old[j * 2 + i] = VGS_old_t[i];
		QCE_old[j * 2 + i] = QCE_old_t[i];
		QCG_old[j * 2 + i] = QCG_old_t[i];
		Qeb_old[j * 2 + i] = Qeb_old_t[i];
	}
	/*
	dVg=Vg-Vg_old[j];//?????????????????? error corrected
	if ((dVg>=10.0)||(dVg<=-10.0))
	{
	transient[j]=1.0;
	}
	submodule(dt,Tvj,Vg,Js,&tcount[j],&ifvi[j],&vrvi[j],&irvi[j],&vfvi[j],&vi[j],&r1[j],&r2[j],&transient[j],
	&Vg_old[j],&ki[j],&kv[j],&Vsw1[j],&Isw1[j],&Vsw2[j],&Isw2[j]);
	Vc[j]=Vsw1[j]+Vsw2[j];
	Vp[j]=Vsw2[j];
	*/
	/*
	if (Vg>=5)
	{
	Gup=1e3;
	Gdn=1e-6;
	}
	else////(Vg<=5)
	{
	Gup=1e-6;
	Gdn=1e3;
	}

	//submodule(Gup,Gdn,&vi[j],Js,&Vc[j],&Vp[j]);

	G11=1/ZC+Gup;
	G12=-Gup;
	G22=Gdn+Gup;
	f1=2*vi[j]/ZC;
	f2=Js;
	Vnew1=(G12*f2-G22*f1)/(G12*G12-G11*G22);
	Vnew2=(G12*f1-G11*f2)/(G12*G12-G11*G22);
	vr=Vnew1-vi[j];
	vi[j]=vr;
	Vc[j]=Vnew1;
	Vp[j]=Vnew2;
	*/
	__syncthreads();

}

__global__ void SM2(int t, double* Vau, double* Uref, double* Js1, double* Js2, double* vi, double* car, double* Vp, double* Vc, double* Vdc,
	double* U, double* Ice1, double* Ice2, double* out1)
{//48D
	int j = threadIdx.x + blockIdx.x * blockDim.x;

	//constants
	double dt = 1e-5;
	//double C=1e-3,ZC=dt/(2*C);//sm
	//double K1=0.5,K2=1,K3=0.5,K4=1;

	double K5, K6, Bu, m, Vg[2];
	double Iu, Id, Js;
	//double vr; 

	//double dVg,Tvj=25.0;
	//simulation starts

	int jcar = j % NoC;
	int phaseshift = 10000 / NoC;
	car[jcar] = 1.0 * abs(((t + 25000 - phaseshift * jcar) % 10000) - 5000) / 5000;

	//-----------------------------------------submodule No. 48
	int jleg = j / NoM;//0,0,0,...,0;1,1,1,...1,1;2,2,2,...,2,2;...;5,5,5,5,5,5,5,5;

	Iu = Js1[jleg];
	Id = Js2[jleg];

	if (Iu >= 0)
		K5 = 0.35 / 20;
	else
		K5 = -0.35 / 20;
	if (Id >= 0)
		K6 = 0.35 / 20;
	else
		K6 = -0.35 / 20;

	int jmmc = j / (3 * NoM);//0,0,0,...,0,0;1,1,1,...,1,1;
	double Vcref = 1.0 * Vdc[jmmc] / NoC;//Vcref=1.0*200.0/NoC;//

	int jvc = j % NoM;//0,1,2,3,4,5,6,7;0,1,2,3,4,5,6,7;...
	if (jvc < NoC)
	{
		Bu = K5 * (Vcref - Vc[j]);
	}
	else
	{
		Bu = K6 * (Vcref - Vc[j]);
	}

	//jleg
	if (jvc < NoC)
	{
		m = (Vau[jleg] + Bu + 0.5 * Vcref - Uref[jleg] / NoC) / Vcref;
	}
	else
	{
		m = (Vau[jleg] + Bu + 0.5 * Vcref + Uref[jleg] / NoC) / Vcref;
	}

	if (jvc < NoC)
	{
		Js = Js1[jleg];
	}
	else
	{
		Js = Js2[jleg];
	}

	if (m > car[jcar] + 0.05)
	{
		Vg[0] = 15;
		Vg[1] = -10;
	}
	else if (m < car[jcar] - 0.05)
	{
		Vg[0] = -10;
		Vg[1] = 15;
	}
	else
	{
		Vg[0] = -10;
		Vg[1] = -10;
	}

	if (j == 3)
	{
		*out1 = Vg[1];
	}

	double U_t[Mb];
	for (int i = 0; i < Mb; i++)
	{
		U_t[i] = U[j * Mb + i];
	}


	submodulesimple(dt, U_t, Vg, &Vc[j], &vi[j], &Vp[j], Js, &Ice1[j], &Ice2[j]);
	for (int i = 0; i < Mb; i++)
	{
		U[j * Mb + i] = U_t[i];
	}


	__syncthreads();

}

__global__ void SM3(int t, double* Uref, double* Js1, double* Js2, double* vi, double* Vp, double* Vc, double* Vdc,
	double* U, double* Ice1, double* Ice2, double* out1, double* Vgg)
{//48D
	int j = threadIdx.x + blockIdx.x * blockDim.x;

	//constants
	double dt = 1e-8;
	//double C=1e-3,ZC=dt/(2*C);//sm
	//double K1=0.5,K2=1,K3=0.5,K4=1;

	double Vg[2];
	double Js;

	Vg[0] = Vgg[2 * j];
	Vg[1] = Vgg[2 * j + 1];

	//-----------------------------------------submodule No. 48
	int jleg = j / NoM;//0,0,0,...,0;1,1,1,...1,1;2,2,2,...,2,2;...;5,5,5,5,5,5,5,5;





	//int jmmc = j / (3 * NoM);//0,0,0,...,0,0;1,1,1,...,1,1;
	//double Vcref = Vdc[jmmc] / NoC;//Vcref=1.0*200.0/NoC;//

	int jvc = j % NoM;//0,1,2,3,4,5,6,7;0,1,2,3,4,5,6,7;...



	if (jvc < NoC)
	{
		Js = Js1[jleg];
	}
	else
	{
		Js = Js2[jleg];
	}








	if (j == 3)
	{
		*out1 = Vg[1];
	}

	double U_t[Mb];
	for (int i = 0; i < Mb; i++)
	{
		U_t[i] = U[j * Mb + i];
	}


	submodulesimple(dt, U_t, Vg, &Vc[j], &vi[j], &Vp[j], Js, &Ice1[j], &Ice2[j]);
	for (int i = 0; i < Mb; i++)
	{
		U[j * Mb + i] = U_t[i];
	}


	__syncthreads();

}

__global__ void SMsumup(double* Vp, double* Vp_up, double* Vp_dn)
{
	int idx = threadIdx.x + blockIdx.x * blockDim.x;//6D, each arm is a dimension

												  //----------------------------------------------------------------sum up
												  /*
												  Vp_up[idx]=0;
												  Vp_dn[idx]=0;
												  for (int i=0;i<NoC;i++)
												  {
												  Vp_up[idx]=Vp_up[idx]+Vp[i+idx*NoM];
												  Vp_dn[idx]=Vp_dn[idx]+Vp[i+NoC+idx*NoM];
												  }
												  */
												  //Vp_up[idx]=Vp[0+idx*NoM]+Vp[1+idx*NoM]+Vp[2+idx*NoM]+Vp[3+idx*NoM];
												  //Vp_dn[idx]=Vp[4+idx*NoM]+Vp[5+idx*NoM]+Vp[6+idx*NoM]+Vp[7+idx*NoM];

	int kk = NoC / 2;
	while (kk != 0)
	{
		for (int k = 0; k < kk; k++)
		{
			Vp[k + idx * NoM] += Vp[k + idx * NoM + kk];
			Vp[k + NoC + idx * NoM] += Vp[k + NoC + idx * NoM + kk];
		}
		__syncthreads();
		kk /= 2;
	}
	Vp_up[idx] = Vp[idx * NoM];
	Vp_dn[idx] = Vp[NoC + idx * NoM];

	__syncthreads();

}

__global__ void SMsumup2(double* Vp, double* Vp_up, double* Vp_dn)
{
	int idx = threadIdx.x + blockIdx.x * blockDim.x;//6D, each arm is a dimension


	Vp_up[idx] = 0;
	Vp_dn[idx] = 0;
	for (int i = 0; i < NoC; i++)
	{
		Vp_up[idx] = Vp_up[idx] + Vp[i + idx * NoM];
		Vp_dn[idx] = Vp_dn[idx] + Vp[i + NoC + idx * NoM];
	}

	//Vp_up[idx]=Vp[0+idx*NoM]+Vp[1+idx*NoM]+Vp[2+idx*NoM]+Vp[3+idx*NoM];
	//Vp_dn[idx]=Vp[4+idx*NoM]+Vp[5+idx*NoM]+Vp[6+idx*NoM]+Vp[7+idx*NoM];											  



	__syncthreads();

}

__global__ void MMCmain(int t, double* theta, double* vsi, double* vsr, double* Vp_up, double* Vp_dn,
	double* Js1, double* Js2, double* vi9, double* vi10, double* vi9pb, double* vi10pb, double* vi9pc, double* vi10pc, double* vLi,
	double* Idaxis, double* Iqaxis, double* vcui, double* vcdi, double* Vdc, double* vilaca, double* vilacb, double* vilacc)//, double *vilaca, double *vilacb, double *vilacc
{//2D
	int idx = threadIdx.x + blockIdx.x * blockDim.x;
	double dt = 1e-8;

	double CC = 10e-3, ZCC = dt / (2 * CC), Zcouple = dt / (1e-6);
	double Vastn, Vbstn, Vcstn, Rarm = 0.0801, Lud = 0.0255, ZLud = 2 * Lud / dt, L = 10e-3, R = 1, Ldc = 1e-6, ZLdc = 2 * Ldc / dt;
	double Lac = 10e-3, ZLac = 2 * Lac / dt;//add

	double Z11, Z12, Z13, Z14, Z15, Z16, U1, Z21, Z22, Z23, Z24, Z25, Z26, U2, Z31, Z32, Z33, Z34, Z35, Z36, U3, Z41, Z42, Z43, Z44, Z45, Z46, U4, Z51, Z52, Z53, Z54, Z55, Z56, U5,
		Z61, Z62, Z63, Z64, Z65, Z66, U6;
	double I[M];
	double Ipa, Ipb, Ipc, Ialpha, Ibeta;

	Vastn = 31e3 / sqrtf(3) * sqrtf(2) * cos(theta[idx]);///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	Vbstn = 31e3 / sqrtf(3) * sqrtf(2) * cos(theta[idx] + 4.0 * PI / 3);///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	Vcstn = 31e3 / sqrtf(3) * sqrtf(2) * cos(theta[idx] + 2.0 * PI / 3);///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
											   //------------------------------------MMC main circuit
	Z11 = 2 * Rarm + 2 * ZCC + 2 * ZLud;
	Z12 = -ZLud - Rarm;
	Z13 = -ZLud - Rarm;
	Z14 = 0;
	Z15 = 0;
	Z16 = -2 * ZCC;
	Z21 = -ZLud - Rarm;
	//Z22 = 2 * ZLud + 2 * Rarm;
	Z22 = 2 * ZLud + 2 * Rarm + 2 * ZLac;//add
	//Z23 = 0;
	Z23 = -2 * ZLac;//add
	//Z24 = -ZLud - Rarm;
	Z24 = -ZLud - Rarm - ZLac;//add
	//Z25 = 0;
	Z25 = ZLac;//add
	Z26 = 0;
	Z31 = -ZLud - Rarm;
	//Z32 = 0;
	Z32 = -2 * ZLac;
	//Z33 = 2 * ZLud + 2 * Rarm;
	Z33 = 2 * ZLud + 2 * Rarm + 2 * ZLac;
	//Z34 = 0;
	Z34 = ZLac;
	//Z35 = -ZLud - Rarm;
	Z35 = -ZLud - Rarm - ZLac;
	Z36 = 0;
	Z41 = 0;
	//Z42 = -ZLud - Rarm;
	Z42 = -ZLud - Rarm - ZLac;
	//Z43 = 0;
	Z43 = ZLac;
	//Z44 = 2 * ZLud + 2 * Rarm;
	Z44 = 2 * ZLud + 2 * Rarm + 2 * ZLac;
	//Z45 = 0;
	Z45 = -2 * ZLac;
	Z46 = 0;
	Z51 = 0;
	//Z52 = 0;
	Z52 = ZLac;
	//Z53 = -ZLud - Rarm;
	Z53 = -ZLud - Rarm - ZLac;
	//Z54 = 0;
	Z54 = -2 * ZLac;
	//Z55 = 2 * ZLud + 2 * Rarm;
	Z55 = 2 * ZLud + 2 * Rarm + 2 * ZLac;
	Z56 = 0;
	Z61 = -2 * ZCC;
	Z62 = 0;
	Z63 = 0;
	Z64 = 0;
	Z65 = 0;
	Z66 = 2 * ZCC + R / 2 + ZLdc + Zcouple;
	U1 = Vp_up[2 + 3 * idx] + Vp_dn[2 + 3 * idx] + 2 * vi9pc[idx] + 2 * vi10pc[idx] - 2 * vcui[idx] - 2 * vcdi[idx];
	//U2 = Vp_up[1 + 3 * idx] - Vp_up[2 + 3 * idx] + 2 * vi9pb[idx] - 2 * vi9pc[idx] + Vbstn - Vcstn;
	U2 = Vp_up[1 + 3 * idx] - Vp_up[2 + 3 * idx] + 2 * vi9pb[idx] - 2 * vi9pc[idx] + Vbstn + 2 * vilacb[idx] - Vcstn - 2 * vilacc[idx];//add
	//U3 = Vp_dn[1 + 3 * idx] - Vp_dn[2 + 3 * idx] + 2 * vi10pb[idx] - 2 * vi10pc[idx] - Vbstn + Vcstn;
	U3 = Vp_dn[1 + 3 * idx] - Vp_dn[2 + 3 * idx] + 2 * vi10pb[idx] - 2 * vi10pc[idx] - Vbstn - 2 * vilacb[idx] + Vcstn + 2 * vilacc[idx];//add
	//U4 = Vp_up[3 * idx] - Vp_up[1 + 3 * idx] + 2 * vi9[idx] - 2 * vi9pb[idx] + Vastn - Vbstn;
	U4 = Vp_up[3 * idx] - Vp_up[1 + 3 * idx] + 2 * vi9[idx] - 2 * vi9pb[idx] + Vastn + 2 * vilaca[idx] - Vbstn - 2 * vilacb[idx];
	//U5 = Vp_dn[3 * idx] - Vp_dn[1 + 3 * idx] + 2 * vi10[idx] - 2 * vi10pb[idx] - Vastn + Vbstn;
	U5 = Vp_dn[3 * idx] - Vp_dn[1 + 3 * idx] + 2 * vi10[idx] - 2 * vi10pb[idx] - Vastn - 2 * vilaca[idx] + Vbstn + 2 * vilacb[idx];
	U6 = 2 * vcui[idx] + 2 * vcdi[idx] - 2 * vLi[idx] - 2 * vsi[idx];

	double AMatrix[M][N] = { Z11,Z12,Z13,Z14,Z15,Z16,U1,Z21,Z22,Z23,Z24,Z25,Z26,U2,Z31,Z32,Z33,Z34,Z35,Z36,U3,Z41,Z42,Z43,Z44,Z45,Z46,U4,Z51,Z52,Z53,Z54,Z55,Z56,U5,Z61,Z62,Z63,Z64,Z65,Z66,U6 };

	solver(AMatrix, I);

	Js1[3 * idx] = -I[3];
	Js2[3 * idx] = -I[4];
	Js1[1 + 3 * idx] = I[3] - I[1];
	Js2[1 + 3 * idx] = I[4] - I[2];
	Js1[2 + 3 * idx] = I[1] - I[0];
	Js2[2 + 3 * idx] = I[2] - I[0];

	double Vnom_sec = 31e3;
	double Pnom = 18e6;

	Ipa = Js1[3 * idx] - Js2[3 * idx];
	Ipb = Js1[1 + 3 * idx] - Js2[1 + 3 * idx];
	Ipc = Js1[2 + 3 * idx] - Js2[2 + 3 * idx];

	Ipa = Ipa * sqrtf(3) * Vnom_sec / sqrtf(2) / Pnom;
	Ipb = Ipb * sqrtf(3) * Vnom_sec / sqrtf(2) / Pnom;
	Ipc = Ipc * sqrtf(3) * Vnom_sec / sqrtf(2) / Pnom;

	Ialpha = (2 * Ipa - Ipb - Ipc) / 3;
	Ibeta = sqrt(3.0) / 3 * (Ipb - Ipc);
	Idaxis[idx] = Ialpha * cosf(theta[idx]) + Ibeta * sinf(theta[idx]);
	Iqaxis[idx] = Ibeta * cosf(theta[idx]) - Ialpha * sinf(theta[idx]);//d,q component in p.u.


	Vdc[idx] = (I[0] - I[5]) * ZCC + 2 * vcui[idx] + (I[0] - I[5]) * ZCC + 2 * vcdi[idx];
	vcui[idx] = (I[0] - I[5]) * ZCC + vcui[idx];
	vcdi[idx] = (I[0] - I[5]) * ZCC + vcdi[idx];
	//the following is newly added for DC yard
	vsr[idx] = I[5] * Zcouple + vsi[idx];
	vLi[idx] = -I[5] * ZLdc - vLi[idx];

	vi9[idx] = -Js1[3 * idx] * ZLud - vi9[idx];
	vi10[idx] = -Js2[3 * idx] * ZLud - vi10[idx];
	vi9pb[idx] = -Js1[1 + 3 * idx] * ZLud - vi9pb[idx];
	vi10pb[idx] = -Js2[1 + 3 * idx] * ZLud - vi10pb[idx];
	vi9pc[idx] = -Js1[2 + 3 * idx] * ZLud - vi9pc[idx];
	vi10pc[idx] = -Js2[2 + 3 * idx] * ZLud - vi10pc[idx];

	vilaca[idx] = -(I[4] - I[3]) * ZLac - vilaca[idx];//add
	vilacb[idx] = -(-I[4] + I[2] + I[3] - I[1]) * ZLac - vilacb[idx];//add
	vilacc[idx] = -(I[1] - I[2]) * ZLac - vilacc[idx];//add
	__syncthreads();
}

__global__ void MMCmain2(int t, double* theta, double* vsi, double* vsr, double* Vp_up, double* Vp_dn,
	double* Js1, double* Js2, double* vi9, double* vi10, double* vi9pb, double* vi10pb, double* vi9pc, double* vi10pc, double* vLi,
	double* Idaxis, double* Iqaxis, double* vcui, double* vcdi, double* Vdc, double* vilaca, double* vilacb, double* vilacc)//, double *vilaca, double *vilacb, double *vilacc
{//2D
	int idx = threadIdx.x + blockIdx.x * blockDim.x;
	double dt = 1e-8;

	double CC = 10e-3, ZCC = dt / (2 * CC), Zcouple = dt / (1e-6);
	double Vastn, Vbstn, Vcstn, Rarm = 0.0801, Lud = 0.0255, ZLud = 2 * Lud / dt, L = 10e-3, R = 1, Ldc = 1e-6, ZLdc = 2 * Ldc / dt;
	double Lac = 10e-3, ZLac = 2 * Lac / dt;//add

	double Z11, Z12, Z13, Z14, Z15, Z16, U1, Z21, Z22, Z23, Z24, Z25, Z26, U2, Z31, Z32, Z33, Z34, Z35, Z36, U3, Z41, Z42, Z43, Z44, Z45, Z46, U4, Z51, Z52, Z53, Z54, Z55, Z56, U5,
		Z61, Z62, Z63, Z64, Z65, Z66, U6;
	double I[M];
	double Ipa, Ipb, Ipc, Ialpha, Ibeta;

	Vastn = 31e3 / sqrtf(3) * sqrtf(2) * cos(theta[idx]);///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	Vbstn = 31e3 / sqrtf(3) * sqrtf(2) * cos(theta[idx] + 4.0 * PI / 3);///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	Vcstn = 31e3 / sqrtf(3) * sqrtf(2) * cos(theta[idx] + 2.0 * PI / 3);///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
											   //------------------------------------MMC main circuit
	Z11 = 2 * Rarm + 2 * ZCC + 2 * ZLud;
	Z12 = -ZLud - Rarm;
	Z13 = -ZLud - Rarm;
	Z14 = 0;
	Z15 = 0;
	Z16 = -2 * ZCC;
	Z21 = -ZLud - Rarm;
	//Z22 = 2 * ZLud + 2 * Rarm;
	Z22 = 2 * ZLud + 2 * Rarm + 2 * ZLac;//add
	//Z23 = 0;
	Z23 = -2 * ZLac;//add
	//Z24 = -ZLud - Rarm;
	Z24 = -ZLud - Rarm - ZLac;//add
	//Z25 = 0;
	Z25 = ZLac;//add
	Z26 = 0;
	Z31 = -ZLud - Rarm;
	//Z32 = 0;
	Z32 = -2 * ZLac;
	//Z33 = 2 * ZLud + 2 * Rarm;
	Z33 = 2 * ZLud + 2 * Rarm + 2 * ZLac;
	//Z34 = 0;
	Z34 = ZLac;
	//Z35 = -ZLud - Rarm;
	Z35 = -ZLud - Rarm - ZLac;
	Z36 = 0;
	Z41 = 0;
	//Z42 = -ZLud - Rarm;
	Z42 = -ZLud - Rarm - ZLac;
	//Z43 = 0;
	Z43 = ZLac;
	//Z44 = 2 * ZLud + 2 * Rarm;
	Z44 = 2 * ZLud + 2 * Rarm + 2 * ZLac;
	//Z45 = 0;
	Z45 = -2 * ZLac;
	Z46 = 0;
	Z51 = 0;
	//Z52 = 0;
	Z52 = ZLac;
	//Z53 = -ZLud - Rarm;
	Z53 = -ZLud - Rarm - ZLac;
	//Z54 = 0;
	Z54 = -2 * ZLac;
	//Z55 = 2 * ZLud + 2 * Rarm;
	Z55 = 2 * ZLud + 2 * Rarm + 2 * ZLac;
	Z56 = 0;
	Z61 = -2 * ZCC;
	Z62 = 0;
	Z63 = 0;
	Z64 = 0;
	Z65 = 0;
	Z66 = 2 * ZCC + R / 2 + ZLdc + Zcouple;
	U1 = Vp_up[2 + 3 * idx] + Vp_dn[2 + 3 * idx] + 2 * vi9pc[idx] + 2 * vi10pc[idx] - 2 * vcui[idx] - 2 * vcdi[idx];
	//U2 = Vp_up[1 + 3 * idx] - Vp_up[2 + 3 * idx] + 2 * vi9pb[idx] - 2 * vi9pc[idx] + Vbstn - Vcstn;
	U2 = Vp_up[1 + 3 * idx] - Vp_up[2 + 3 * idx] + 2 * vi9pb[idx] - 2 * vi9pc[idx] + Vbstn + 2 * vilacb[idx] - Vcstn - 2 * vilacc[idx];//add
	//U3 = Vp_dn[1 + 3 * idx] - Vp_dn[2 + 3 * idx] + 2 * vi10pb[idx] - 2 * vi10pc[idx] - Vbstn + Vcstn;
	U3 = Vp_dn[1 + 3 * idx] - Vp_dn[2 + 3 * idx] + 2 * vi10pb[idx] - 2 * vi10pc[idx] - Vbstn - 2 * vilacb[idx] + Vcstn + 2 * vilacc[idx];//add
	//U4 = Vp_up[3 * idx] - Vp_up[1 + 3 * idx] + 2 * vi9[idx] - 2 * vi9pb[idx] + Vastn - Vbstn;
	U4 = Vp_up[3 * idx] - Vp_up[1 + 3 * idx] + 2 * vi9[idx] - 2 * vi9pb[idx] + Vastn + 2 * vilaca[idx] - Vbstn - 2 * vilacb[idx];
	//U5 = Vp_dn[3 * idx] - Vp_dn[1 + 3 * idx] + 2 * vi10[idx] - 2 * vi10pb[idx] - Vastn + Vbstn;
	U5 = Vp_dn[3 * idx] - Vp_dn[1 + 3 * idx] + 2 * vi10[idx] - 2 * vi10pb[idx] - Vastn - 2 * vilaca[idx] + Vbstn + 2 * vilacb[idx];
	U6 = 2 * vcui[idx] + 2 * vcdi[idx] - 2 * vLi[idx] - 2 * vsi[idx];

	double AMatrix[M][N] = { Z11,Z12,Z13,Z14,Z15,Z16,U1,Z21,Z22,Z23,Z24,Z25,Z26,U2,Z31,Z32,Z33,Z34,Z35,Z36,U3,Z41,Z42,Z43,Z44,Z45,Z46,U4,Z51,Z52,Z53,Z54,Z55,Z56,U5,Z61,Z62,Z63,Z64,Z65,Z66,U6 };

	solver(AMatrix, I);

	Js1[3 * idx] = -I[3];
	Js2[3 * idx] = -I[4];
	Js1[1 + 3 * idx] = I[3] - I[1];
	Js2[1 + 3 * idx] = I[4] - I[2];
	Js1[2 + 3 * idx] = I[1] - I[0];
	Js2[2 + 3 * idx] = I[2] - I[0];

	double Vnom_sec = 31e3;
	double Pnom = 18e6;

	Ipa = Js1[3 * idx] - Js2[3 * idx];
	Ipb = Js1[1 + 3 * idx] - Js2[1 + 3 * idx];
	Ipc = Js1[2 + 3 * idx] - Js2[2 + 3 * idx];

	Ipa = Ipa * sqrtf(3) * Vnom_sec / sqrtf(2) / Pnom;
	Ipb = Ipb * sqrtf(3) * Vnom_sec / sqrtf(2) / Pnom;
	Ipc = Ipc * sqrtf(3) * Vnom_sec / sqrtf(2) / Pnom;

	Ialpha = (2 * Ipa - Ipb - Ipc) / 3;
	Ibeta = sqrt(3.0) / 3 * (Ipb - Ipc);
	Idaxis[idx] = Ialpha * cosf(theta[idx]) + Ibeta * sinf(theta[idx]);
	Iqaxis[idx] = Ibeta * cosf(theta[idx]) - Ialpha * sinf(theta[idx]);


	Vdc[idx] = (I[0] - I[5]) * ZCC + 2 * vcui[idx] + (I[0] - I[5]) * ZCC + 2 * vcdi[idx];
	vcui[idx] = (I[0] - I[5]) * ZCC + vcui[idx];
	vcdi[idx] = (I[0] - I[5]) * ZCC + vcdi[idx];
	//the following is newly added for DC yard
	vsr[idx] = I[5] * Zcouple + vsi[idx];
	vLi[idx] = -I[5] * ZLdc - vLi[idx];

	vi9[idx] = -Js1[3 * idx] * ZLud - vi9[idx];
	vi10[idx] = -Js2[3 * idx] * ZLud - vi10[idx];
	vi9pb[idx] = -Js1[1 + 3 * idx] * ZLud - vi9pb[idx];
	vi10pb[idx] = -Js2[1 + 3 * idx] * ZLud - vi10pb[idx];
	vi9pc[idx] = -Js1[2 + 3 * idx] * ZLud - vi9pc[idx];
	vi10pc[idx] = -Js2[2 + 3 * idx] * ZLud - vi10pc[idx];

	vilaca[idx] = -(I[4] - I[3]) * ZLac - vilaca[idx];//add
	vilacb[idx] = -(-I[4] + I[2] + I[3] - I[1]) * ZLac - vilacb[idx];//add
	vilacc[idx] = -(I[1] - I[2]) * ZLac - vilacc[idx];//add
	__syncthreads();
}


__device__ void datasort2(double Vxy[NoSTN], double Ixy[NoSTN], double Ik[NLine], double Im[NLine], double Vk[NLine], double Vm[NLine])
{
	//int idx = threadIdx.x+blockIdx.x*blockDim.x;//1 thread

	Vk[0] = Vxy[0];//left: line #; right: MMC #


	Vm[0] = Vxy[1];//left: line #; right: MMC #


	Ik[0] = Ixy[0];//left: line #; right: MMC branch #


	Im[0] = Ixy[1];//left: line #; right: MMC branch #


	__syncthreads();
}


__device__ void datasort1(double* Ihisk, double* Ihism, double* Ihisxy)
{
	//int idx = threadIdx.x+blockIdx.x*blockDim.x;//1 thread

	Ihisxy[0] = Ihisk[0];//left:MMC0, right:L0
	Ihisxy[1] = Ihism[0];//left:MMC0, right:L1

	__syncthreads();
}

//FDLM
__global__ void BurgeronLM(double* Vxy, double* Ixy, double* Ihisxy)//(double *Ik,double *Im,double *Vk,double *Vm,double *Ihisk,double *Ihism)
{
	int idx = threadIdx.x + blockIdx.x * blockDim.x;//1*NLine threads

	double Rdis = 0.0121 * 50, Ldis = 0.1056e-3, Cdis = 0.2961e-6, Z0 = sqrtf(Ldis / Cdis), Z00 = Z0 + Rdis / 4, H = (Z0 - Rdis / 4) / (Z0 + Rdis / 4);

	__shared__ double Ik[NLine], Im[NLine], Vk[NLine], Vm[NLine];
	if (idx == 0)
	{
		datasort2(Vxy, Ixy, Ik, Im, Vk, Vm);
	}

	__shared__ double Ihisk[NLine], Ihism[NLine];
	Ihisk[idx] = Vm[idx] / Z00 + H * Im[idx];
	Ihism[idx] = Vk[idx] / Z00 + H * Ik[idx];


	if (idx == 0)
	{
		datasort1(Ihisk, Ihism, Ihisxy);
	}

	__syncthreads();
}




//MTDC system changes made

__global__ void DCyard(double* vmi, double* vmr, double* vixy, double* Jsxy, double* Ihisxy, double* Vxy, double* Ixy, double* obs)
{
	int jj = threadIdx.x;//11 DCyards corresponds to 11 MMCs; 0,1,2,3,4,5,6,7,8,9,10

	double dt = 1e-8;


	double Rdis = 0.0121 * 50, Ldis = 0.1056e-3, Cdis = 0.2961e-6, Z0 = sqrtf(Ldis / Cdis), Z00 = Z0 + Rdis / 4, H = (Z0 - Rdis / 4) / (Z0 + Rdis / 4);

	double L = 100e-3, ZL = 2 * L / dt, ZTL = Z0 + Rdis / 4, Zcouple = dt / (1e-6), r = 0.5;//Zcouple for DC yard coupling, ZTL is for transmission line coupling


	double Z11, Z12, Z22, U1, U2, Im1, Im2;

	Z11 = Zcouple + r + ZL + ZTL;
	Z12 = 0;
	Z22 = 0;
	U1 = 2 * vmi[jj] - 2 * vixy[jj] - Ihisxy[jj] * Z00;
	U2 = 0;
	Im1 = U1 / Z11;
	Im2 = 0;
	vmr[jj] = -Im1 * Zcouple + vmi[jj];


	vixy[jj] = -Im1 * ZL - vixy[jj];
	Jsxy[jj] = Im1;


	Vxy[jj] = Im1 * Z00 + Ihisxy[jj] * Z00;

	Ixy[jj] = Vxy[jj] / Z00 - Ihisxy[jj];

	if (jj == 0)
	{
		*obs = Im1;
	}
}

__global__ void DCyard2(double* vmi, double* vmr, double* vixy, double* Jsxy, double* Ihisxy, double* Vxy, double* Ixy, double* obs)
{
	int jj = threadIdx.x;//11 DCyards corresponds to 11 MMCs; 0,1,2,3,4,5,6,7,8,9,10

	double dt = 1e-8;


	double Rdis = 0.0121 * 50, Ldis = 0.1056e-3, Cdis = 0.2961e-6, Z0 = sqrtf(Ldis / Cdis), Z00 = Z0 + Rdis / 4, H = (Z0 - Rdis / 4) / (Z0 + Rdis / 4);

	double L = 100e-3, ZL = 2 * L / dt, ZTL = Z0 + Rdis / 4, Zcouple = dt / (1e-6), r = 0.5;//Zcouple for DC yard coupling, ZTL is for transmission line coupling


	double Z11, Z12, Z22, U1, U2, Im1, Im2;

	Z11 = Zcouple + r + ZL + ZTL;
	Z12 = 0;
	Z22 = 0;
	U1 = 2 * vmi[jj] - 2 * vixy[jj] - Ihisxy[jj] * Z00;
	U2 = 0;
	Im1 = U1 / Z11;
	Im2 = 0;
	vmr[jj] = -Im1 * Zcouple + vmi[jj];


	vixy[jj] = -Im1 * ZL - vixy[jj];
	Jsxy[jj] = Im1;


	Vxy[jj] = Im1 * Z00 + Ihisxy[jj] * Z00;

	Ixy[jj] = Vxy[jj] / Z00 - Ihisxy[jj];

	if (jj == 0)
	{
		*obs = Im1;
	}
}

__global__ void HVDC(double* output1, double* output2, double* output3, double* output4, double* output5, double* output6)//
{
	int idex = threadIdx.x;//+blockIdx.x*blockDim.x;
	double dt = 1e-8;

	double* Vdc = { 0 }, * Js1 = { 0 }, * Js2 = { 0 }, * Vzu2 = { 0 }, * Vau2 = { 0 }, * Vau = { 0 };
	cudaMalloc((void**)&Vdc, NoSTN * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&Js1, 3 * NoSTN * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&Js2, 3 * NoSTN * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&Vzu2, 3 * NoSTN * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&Vau2, 3 * NoSTN * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&Vau, 3 * NoSTN * sizeof(double));//don's forget memory allocation

	double* vi9 = { 0 }, * vi10 = { 0 }, * vi9pb = { 0 }, * vi10pb = { 0 }, * vi9pc = { 0 }, * vi10pc = { 0 }, * vLi = { 0 };
	double* vilaca = { 0 }, * vilacb = { 0 }, * vilacc = { 0 };//add
	cudaMalloc((void**)&vi9, NoSTN * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&vi10, NoSTN * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&vi9pb, NoSTN * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&vi10pb, NoSTN * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&vi9pc, NoSTN * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&vi10pc, NoSTN * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&vLi, NoSTN * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&vilaca, NoSTN * sizeof(double));//add
	cudaMalloc((void**)&vilacb, NoSTN * sizeof(double));//add
	cudaMalloc((void**)&vilacc, NoSTN * sizeof(double));//add

	double* theta = { 0 }, * Uref = { 0 }, * sum1 = { 0 }, * sum2 = { 0 }, * Idaxis = { 0 }, * Iqaxis = { 0 }, * sum3 = { 0 }, * sum4 = { 0 }, * sum5 = { 0 }, * vcui = { 0 }, * vcdi = { 0 };
	cudaMalloc((void**)&theta, NoSTN * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&Uref, 3 * NoSTN * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&sum1, NoSTN * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&sum2, NoSTN * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&Idaxis, NoSTN * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&Iqaxis, NoSTN * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&sum3, NoSTN * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&sum4, NoSTN * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&sum5, NoSTN * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&vcui, NoSTN * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&vcdi, NoSTN * sizeof(double));//don's forget memory allocation
	/////////////////////////////////////////////////////////////////////////////////////////
	for (int i = 0; i < NoSTN; i++)
	{
		vcui[i] = 15e3;
		vcdi[i] = 15e3;
	}
	////////////////////////////////////////////////////////////////////////////////////////////
													   //
	double* Vp_up = { 0 }, * Vp_dn = { 0 }, * vi = { 0 }, * Vp = { 0 }, * Vc = { 0 };
	cudaMalloc((void**)&Vp_up, 3 * NoSTN * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&Vp_dn, 3 * NoSTN * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&vi, 3 * NoSTN * NoM * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&Vp, 3 * NoSTN * NoM * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&Vc, 3 * NoSTN * NoM * sizeof(double));//don's forget memory allocation

	double* car = { 0 };
	cudaMalloc((void**)&car, NoC * sizeof(double));//don's forget memory allocation
													//

	double* vsi = { 0 }, * vsr = { 0 }, * vmi = { 0 }, * vmr = { 0 };
	cudaMalloc((void**)&vsi, NoSTN * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&vsr, NoSTN * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&vmi, NoSTN * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&vmr, NoSTN * sizeof(double));//don's forget memory allocation

													  //double *obs1,*obs2;
													  //cudaMalloc((void **) &obs1, sizeof(double));//don's forget memory allocation
													  //cudaMalloc((void **) &obs2, sizeof(double));//don's forget memory allocation


													  //--------------------------------------------------------------------------------HHB
	double* Vhhb = { 0 }, * vixy = { 0 }, * Jhhb = { 0 };
	double* U1 = { 0 }, * U2 = { 0 }, * LCS, * UFD, * MCB, * vi2 = { 0 }, * vi3 = { 0 };
	double* tcount = { 0 }, * trip = { 0 }, * obs;

	cudaMalloc((void**)&Vhhb, NoSTN * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&vixy, NoSTN * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&Jhhb, NoSTN * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&U1, NoSTN * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&U2, NoSTN * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&LCS, NoSTN * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&UFD, NoSTN * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&MCB, NoSTN * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&vi2, NoSTN * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&vi3, NoSTN * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&tcount, NoSTN * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&trip, NoSTN * sizeof(double));//don's forget memory allocation

	cudaMalloc((void**)&obs, sizeof(double));//don's forget memory allocation

	for (int i = 0; i < NoSTN; i++)
	{
		LCS[i] = 1e-3, UFD[i] = 1e-3, MCB[i] = 1e-3;
	}
	//----------------------------------------------------------------Burgeron line model


	double* Ihisxy, * Vxy, * Ixy;/////////////////////////////////////////////////////////
	double* Vk, * Vm, * Ik, * Im, * Ihisk, * Ihism;// *Iki, *Imi, *Ikr, *Imr;
	cudaMalloc((void**)&Vxy, NoSTN * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&Ixy, NoSTN * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&Ihisxy, NoSTN * sizeof(double));//don's forget memory allocation

	cudaMalloc((void**)&Vk, NLine * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&Vm, NLine * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&Ik, NLine * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&Im, NLine * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&Ihisk, NLine * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&Ihism, NLine * sizeof(double));//don's forget memory allocation

	for (int i = 0; i < NoSTN; i++)
	{
		Ihisxy[i] = 0;
	}

	for (int i = 0; i < NLine; i++)
	{
		Ihisk[i] = 0;
		Ihism[i] = 0;
	}



	//////////////////////////////////////////////////////////////////////////////////
	/////////////////////////////////////////////////////////////////////////////////

	double* U_IGBT = { 0 }, * U = { 0 }, * VGS_old, * QCE_old, * QCG_old, * Qeb_old, * Ice1, * Ice2;//Submodule
	double* out1;
	cudaMalloc((void**)&U_IGBT, Ma * 3 * NoSTN * NoM * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&U, Mb * 3 * NoSTN * NoM * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&VGS_old, 2 * 3 * NoSTN * NoM * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&QCE_old, 2 * 3 * NoSTN * NoM * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&QCG_old, 2 * 3 * NoSTN * NoM * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&Qeb_old, 2 * 3 * NoSTN * NoM * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&Ice1, 3 * NoSTN * NoM * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&Ice2, 3 * NoSTN * NoM * sizeof(double));//don's forget memory allocation

	cudaMalloc((void**)&out1, 3 * NoSTN * NoM * sizeof(double));//don's forget memory allocation

	for (int i = 0; i < 3 * NoSTN * NoM; i++)
	{
		vi[i] = 625;
	}
	for (int i = 0; i < 2 * 3 * NoM * NoSTN; i++)
	{
		VGS_old[i] = 0; QCE_old[i] = 0; QCG_old[i] = 0; Qeb_old[i] = 0;
	}
	//for (int i = 0; i<Ma * 3 * NoSTN*NoM; i++)
	//{
	//	U[i] = 0;
	//}
	double* Vg, * Deadtime_count;
	cudaMalloc((void**)&Vg, 2 * 3 * NoM * NoSTN * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&Deadtime_count, 2 * 3 * NoM * NoSTN * sizeof(double));//don's forget memory allocation
	double* NU, * ND;
	cudaMalloc((void**)&NU, 2 * sizeof(double));//don's forget memory allocation
	cudaMalloc((void**)&ND, 2 * sizeof(double));//don's forget memory allocation
	//-----------------------------------------------------------------------------

	int t1 = 0;


	while (t1 < RunTime1)
	{

		PQcontrol << <1, NoSTN >> > (theta, Uref, sum1, sum2, Idaxis, Iqaxis, sum3, sum4, sum5, Vdc);

		//PSC << <1, 3 * NoSTN >> >(Vc, Vzu2, Vau2, Vau, Vdc, Js1, Js2);

		NLM2 << <1, 3 * NoSTN >> > (t1, Vc, Uref, Js1, Js2, Vg, Deadtime_count, NU, ND);

		SM << <56, 3 * NoSTN * NoM / 56 >> > (t1, Uref, Js1, Js2, vi, Vp, Vc, Vdc,
			U_IGBT, VGS_old, QCE_old, QCG_old, Qeb_old, Ice1, Ice2, out1, Vg);
		//SM << <4, 3 * NoSTN*NoM / 4 >> >(t, Vau, Uref, Js1, Js2, vi, car, Vp, Vc, Vdc,
			//U, qcg_old, vcge_old, Ice1, Ice2, out1);//increase the number of blocks here to enable higher MMC level

		//SM << <4, 3 * NoSTN*NoM / 4 >> >(t, Vau, Uref, Js1, Js2, vi, car, Vp, Vc, Vdc,
		//	U, VGS_old, QCE_old, QCG_old, Qeb_old, Ice1, Ice2, out1);//increase the number of blocks here to enable higher MMC level

		//SM2 << <200, 3 * NoSTN*NoM / 200 >> > (t, Vau, Uref, Js1, Js2, vi, car, Vp, Vc, Vdc,
		//	U, Ice1, Ice2, out1);

		SM3 << <56, 3 * NoSTN * NoM / 56 >> > (t1, Uref, Js1, Js2, vi, Vp, Vc, Vdc, U, Ice1, Ice2, out1, Vg);

		SMsumup2 << <1, 3 * NoSTN >> > (Vp, Vp_up, Vp_dn);

		MMCmain << <1, NoSTN >> > (t1, theta, vsi, vsr, Vp_up, Vp_dn, Js1, Js2, vi9, vi10, vi9pb, vi10pb, vi9pc, vi10pc, vLi, Idaxis, Iqaxis, vcui, vcdi, Vdc, vilaca, vilacb, vilacc);//


		DCyard << <1, NoSTN >> > (vsr, vsi, vixy, Jhhb, Ihisxy, Vxy, Ixy, obs);//Js -> Jhhb, Vp -> Vhhb,

		BurgeronLM << <1, NLine >> > (Vxy, Ixy, Ihisxy);

		cudaDeviceSynchronize();


		output1[t1] = Vp[0];
		output2[t1] = Vp[1];
		output3[t1] = (Vp_dn[0] - Vp_up[0]) / 2;

		output4[t1] = Vdc[0];
		output5[t1] = Vdc[1];
		output6[t1] = Vc[0];

		t1 = t1 + 1;
		__syncthreads();
	}
	
	/*int t2 = 0;
	while (t2 < RunTime2)
	{

		PQcontrol2 << <1, NoSTN >> > (theta, Uref, sum1, sum2, Idaxis, Iqaxis, sum3, sum4, sum5, Vdc);

		//PSC << <1, 3 * NoSTN >> >(Vc, Vzu2, Vau2, Vau, Vdc, Js1, Js2);

		NLM22 << <1, 3 * NoSTN >> > (t2, Vc, Uref, Js1, Js2, Vg, Deadtime_count, NU, ND);

		//SM << <4, 3 * NoSTN*NoM / 4 >> >(t, Vau, Uref, Js1, Js2, vi, car, Vp, Vc, Vdc,
			//U, qcg_old, vcge_old, Ice1, Ice2, out1);//increase the number of blocks here to enable higher MMC level

		SM << <56, 3 * NoSTN * NoM / 56 >> > (t2, Uref, Js1, Js2, vi, Vp, Vc, Vdc,
			U_IGBT, VGS_old, QCE_old, QCG_old, Qeb_old, Ice1, Ice2, out1, Vg);//increase the number of blocks here to enable higher MMC level

		//SM2 << <200, 3 * NoSTN*NoM / 200 >> > (t, Vau, Uref, Js1, Js2, vi, car, Vp, Vc, Vdc,
		//	U, Ice1, Ice2, out1);

		//SM3 << <56, 3 * NoSTN*NoM / 56 >> > (t2, Uref, Js1, Js2, vi, Vp, Vc, Vdc, U, Ice1, Ice2, out1, Vg);

		SMsumup2 << <1, 3 * NoSTN >> > (Vp, Vp_up, Vp_dn);

		MMCmain2 << <1, NoSTN >> > (t2, theta, vsi, vsr, Vp_up, Vp_dn, Js1, Js2, vi9, vi10, vi9pb, vi10pb, vi9pc, vi10pc, vLi, Idaxis, Iqaxis, vcui, vcdi, Vdc, vilaca, vilacb, vilacc);//


		DCyard2 << <1, NoSTN >> > (vsr, vsi, vixy, Jhhb, Ihisxy, Vxy, Ixy, obs);//Js -> Jhhb, Vp -> Vhhb,

		BurgeronLM << <1, NLine >> > (Vxy, Ixy, Ihisxy);

		cudaDeviceSynchronize();

		output1[t2] = Vp[0];
		output2[t2] = Vp[1];
		output3[t2] = (Vp_dn[0] - Vp_up[0]) / 2;

		output4[t2] = Vdc[0];
		output5[t2] = Vdc[1];
		output6[t2] = Vc[0];

		//output1[t] = Vp_up[0];
		//output2[t] = Vp_dn[0];
		//output3[t] = 2 * vi9[0] + Js1[0] * (0.0801 + 2 * 0.0255 / 1e-5);

		//output4[t] = 2 * vi10[0] + Js2[0] * (0.0801 + 2 * 0.0255 / 1e-5);
		//output5[t] = Vp_up[1] + 2 * vi9pb[0] + Vp_dn[1] + 2 * vi10pb[0] + Js1[1] * (0.0801 + 2 * 0.0255 / 1e-5) + Js2[1] * (0.0801 + 2 * 0.0255 / 1e-5);
		//output6[t] = Vp_up[0] + 2 * vi9[0] + Vp_dn[0] + 2 * vi10[0] + Js1[0] * (0.0801 + 2 * 0.0255 / 1e-5) + Js2[0] * (0.0801 + 2 * 0.0255 / 1e-5);

		t2 = t2 + 1;
		
		__syncthreads();
	}
	*/
	
}

int main()
{
	FILE* fp1 = fopen("output1.txt", "w");//??????
	FILE* fp2 = fopen("output2.txt", "w");//??????
	FILE* fp3 = fopen("output3.txt", "w");//??????
	if (fp2 == NULL || fp1 == NULL || fp3 == NULL) {//??????????
		puts("!");
		return 0;
	}
	FILE* fp4 = fopen("output4.txt", "w");//??????
	FILE* fp5 = fopen("output5.txt", "w");//??????
	FILE* fp6 = fopen("output6.txt", "w");//??????
	if (fp4 == NULL || fp5 == NULL || fp6 == NULL) {//??????????
		puts("!");
		return 0;
	}

	double* host_out1, * host_out2, * host_out3;
	double* dev_out1, * dev_out2, * dev_out3;
	host_out1 = (double*)malloc(RunTime1 * sizeof(double));
	host_out2 = (double*)malloc(RunTime1 * sizeof(double));
	host_out3 = (double*)malloc(RunTime1 * sizeof(double));
	cudaMalloc((void**)&dev_out1, RunTime1 * sizeof(double));
	cudaMalloc((void**)&dev_out2, RunTime1 * sizeof(double));
	cudaMalloc((void**)&dev_out3, RunTime1 * sizeof(double));
	cudaMemcpy(dev_out1, host_out1, RunTime1 * sizeof(double), cudaMemcpyHostToDevice);
	cudaMemcpy(dev_out2, host_out2, RunTime1 * sizeof(double), cudaMemcpyHostToDevice);
	cudaMemcpy(dev_out3, host_out3, RunTime1 * sizeof(double), cudaMemcpyHostToDevice);

	double* host_out4, * host_out5, * host_out6;
	double* dev_out4, * dev_out5, * dev_out6;
	host_out4 = (double*)malloc(RunTime1 * sizeof(double));
	host_out5 = (double*)malloc(RunTime1 * sizeof(double));
	host_out6 = (double*)malloc(RunTime1 * sizeof(double));
	cudaMalloc((void**)&dev_out4, RunTime1 * sizeof(double));
	cudaMalloc((void**)&dev_out5, RunTime1 * sizeof(double));
	cudaMalloc((void**)&dev_out6, RunTime1 * sizeof(double));
	cudaMemcpy(dev_out4, host_out4, RunTime1 * sizeof(double), cudaMemcpyHostToDevice);
	cudaMemcpy(dev_out5, host_out5, RunTime1 * sizeof(double), cudaMemcpyHostToDevice);
	cudaMemcpy(dev_out6, host_out6, RunTime1 * sizeof(double), cudaMemcpyHostToDevice);

	cudaEvent_t start, stop;//unit: ms

	cudaEventCreate(&start);
	cudaEventCreate(&stop);
	cudaEventRecord(start, 0);

	//GPU
	HVDC <<<1, 1 >>> (dev_out1, dev_out2, dev_out3, dev_out4, dev_out5, dev_out6);//

	cudaEventRecord(stop, 0);
	cudaEventSynchronize(stop);
	float elapsedTime;
	cudaEventElapsedTime(&elapsedTime, start, stop);

	cudaMemcpy(host_out1, dev_out1, RunTime1 * sizeof(double), cudaMemcpyDeviceToHost);
	cudaMemcpy(host_out2, dev_out2, RunTime1 * sizeof(double), cudaMemcpyDeviceToHost);
	cudaMemcpy(host_out3, dev_out3, RunTime1 * sizeof(double), cudaMemcpyDeviceToHost);
	cudaMemcpy(host_out4, dev_out4, RunTime1 * sizeof(double), cudaMemcpyDeviceToHost);
	cudaMemcpy(host_out5, dev_out5, RunTime1 * sizeof(double), cudaMemcpyDeviceToHost);
	cudaMemcpy(host_out6, dev_out6, RunTime1 * sizeof(double), cudaMemcpyDeviceToHost);

	for (int k = 0; k < RunTime1; k++)
	{
		//if (k % 100 == 0)
		//{
		fprintf(fp1, "%d %f\n", k, host_out1[k]);//???????????
		fprintf(fp2, "%d %f\n", k, host_out2[k]);//???????????
		fprintf(fp3, "%d %f\n", k, host_out3[k]);//???????????
		fprintf(fp4, "%d %f\n", k, host_out4[k]);//???????????
		fprintf(fp5, "%d %f\n", k, host_out5[k]);//???????????
		fprintf(fp6, "%d %f\n", k, host_out6[k]);//???????????
												 //printf("%d, Vc = {%f,%f}\n", k, host_O[k], host_O1[k]);
	//}

	}
	printf("%f ms\n", elapsedTime);
	fclose(fp1);//??????,?????
	fclose(fp2);//??????,?????
	fclose(fp3);//??????,?????
	fclose(fp4);//??????,?????
	fclose(fp5);//??????,?????
	fclose(fp6);//??????,?????

	cudaFree(dev_out1);
	cudaFree(dev_out2);
	cudaFree(dev_out3);
	cudaFree(dev_out4);
	cudaFree(dev_out5);
	cudaFree(dev_out6);
	return 0;
}